Product Reviews Clustering

1. Project Introduction

In the fast-goods industry, customer feedback is crucial for maintaining product quality and brand reputation. With thousands of reviews submitted daily on platforms like Sephora, manually analyzing this unstructured text data is impossible.

The motivation and goal for this analysis is to apply unsupervised learning techniques to automatically uncover distinct themes within customer reviews, transforming raw text information into actionable business insights.

The intended audience for this analysis is the product development or other related teams at a beauty company. By identifying specific clusters of customer feedback, the team can take targeted next steps: for instance, if a cluster reveals consistent complaints about “leaking bottles”, the team can redesign the packaging; if another cluster highlights “skin irritation”, the formula can be revisited for safety testing.

The dataset we will be using is Sephora Products and Skincare Reviews from Kaggle, which was collected via Python scraper in March 2023 and contains:

• Product Dataset: information about all beauty products (over 8,000) from the Sephora online store, including product and brand names, prices, ingredients, ratings, and all features.

• Customer Review Dataset: user reviews (about 1 million on over 2,000 products) of all products from the Skincare category, including user appearances, and review ratings by other users.

Given the limited computational resources available and practical purposes for this project, we adopted a two-stage data approach. The complete dataset was used for data cleaning and EDA to identify global trends. For the computationally intensive clustering algorithms, we utilized a subset of approximately 10,000 stratified sample reviews to ensure efficient model training and parameter tuning.

# import relevant libraries

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px


from collections import defaultdict
from wordcloud import STOPWORDS
import string

import ast
import re

import operator
import unicodedata
from collections import Counter

from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split
from umap import UMAP
from sklearn.metrics import calinski_harabasz_score, davies_bouldin_score, silhouette_score
import hdbscan
from sklearn.manifold import TSNE
from sklearn.preprocessing import normalize
from sklearn.cluster import KMeans, MiniBatchKMeans
from sklearn.manifold import trustworthiness
from scipy.spatial.distance import pdist, squareform
from sklearn.neighbors import NearestNeighbors
from sklearn.cluster import DBSCAN
from sklearn.metrics import adjusted_rand_score, silhouette_samples

import warnings
warnings.filterwarnings("ignore", message="MiniBatchKMeans is known to have a memory leak on Windows")

from wordcloud import WordCloud
import math

plt.style.use("seaborn-v0_8-white")

import plotly.io as pio
pio.renderers.default = 'notebook'

2. Data Loading

Load product dataset.

product_df = pd.read_csv("dataset/Sephora/product_info.csv")
product_df.head()

	product_id	product_name	brand_id	brand_name	loves_count	rating	reviews	size	variation_type	variation_value	...	online_only	out_of_stock	sephora_exclusive	highlights	primary_category	secondary_category	tertiary_category	child_count	child_max_price	child_min_price
0	P473671	Fragrance Discovery Set	6342	19-69	6320	3.6364	11.0	NaN	NaN	NaN	...	1	0	0	['Unisex/ Genderless Scent', 'Warm &Spicy Scen...	Fragrance	Value & Gift Sets	Perfume Gift Sets	0	NaN	NaN
1	P473668	La Habana Eau de Parfum	6342	19-69	3827	4.1538	13.0	3.4 oz/ 100 mL	Size + Concentration + Formulation	3.4 oz/ 100 mL	...	1	0	0	['Unisex/ Genderless Scent', 'Layerable Scent'...	Fragrance	Women	Perfume	2	85.0	30.0
2	P473662	Rainbow Bar Eau de Parfum	6342	19-69	3253	4.2500	16.0	3.4 oz/ 100 mL	Size + Concentration + Formulation	3.4 oz/ 100 mL	...	1	0	0	['Unisex/ Genderless Scent', 'Layerable Scent'...	Fragrance	Women	Perfume	2	75.0	30.0
3	P473660	Kasbah Eau de Parfum	6342	19-69	3018	4.4762	21.0	3.4 oz/ 100 mL	Size + Concentration + Formulation	3.4 oz/ 100 mL	...	1	0	0	['Unisex/ Genderless Scent', 'Layerable Scent'...	Fragrance	Women	Perfume	2	75.0	30.0
4	P473658	Purple Haze Eau de Parfum	6342	19-69	2691	3.2308	13.0	3.4 oz/ 100 mL	Size + Concentration + Formulation	3.4 oz/ 100 mL	...	1	0	0	['Unisex/ Genderless Scent', 'Layerable Scent'...	Fragrance	Women	Perfume	2	75.0	30.0

5 rows × 27 columns

2.1 Column Type Inconsistency (Review Dataset)

Load customer review dataset.

review_df_1 = pd.read_csv("dataset/Sephora/reviews_0-250.csv", index_col=0)

C:\Users\fuyuz\AppData\Local\Temp\ipykernel_34616\3989447197.py:1: DtypeWarning:

Columns (1) have mixed types. Specify dtype option on import or set low_memory=False.

By loading the first partial customer review dataset, we see that there are columns having data type inconsistency. We need to check and correct this issue for each partial review dataset.

The problematic columns are author_id and user profile columns, and the primary problem is that some author_ids are stored as string type instead of an integer, every new encountered string author_id will be replaced with an integer, starting from 1. User profile columns such as skin_tone and eye_color contain both NA values and strings. All NA values in these columns have been replaced with ‘No_profile’. Same to the review_text and review_title, all NA values will be replaced with ‘No_review’ and ‘No_review_title’ respectively.

# check data tpye in each column 
def column_type_check(df): # data frame
    print('Columns that have multiple data types: ')

    for column in df.columns:
        if len(df[column].apply(type).value_counts()) >= 2:
            print(' ', column)

2.1.1 First Partial Review Dataset

review_df_1 = pd.read_csv("dataset/Sephora/reviews_0-250.csv", index_col=0)
column_type_check(review_df_1)

C:\Users\fuyuz\AppData\Local\Temp\ipykernel_34616\467310262.py:1: DtypeWarning:

Columns (1) have mixed types. Specify dtype option on import or set low_memory=False.

Columns that have multiple data types: 
  author_id
  review_text
  review_title
  skin_tone
  eye_color
  skin_type
  hair_color

#author_id
mask = (review_df_1.loc[:,'author_id'] == 'dummyUser')
review_df_1.loc[mask,'author_id'] = 1

mask = review_df_1.iloc[:,0].apply(lambda x: isinstance(x, str)) 
a = review_df_1[mask]
mask_idx = a[a.loc[:,'author_id'].str.contains(r'order', regex=False)].index # for the str 'order...'

unique_id = review_df_1.loc[mask_idx, 'author_id'].unique()

# create id mapping for storing id generated 
author_id_mapping = {old: new_id for new_id, old in enumerate(unique_id, start= 2)}

review_df_1.loc[mask_idx, 'author_id'] = review_df_1.loc[mask_idx, 'author_id'].map(author_id_mapping).astype(int)

review_df_1['author_id'] = review_df_1['author_id'].astype(int)

#review_text
mask = review_df_1.iloc[:,8].apply(lambda x: isinstance(x, float))
review_df_1.loc[mask,'review_text'] = 'No_review'

#review_title
mask = review_df_1.iloc[:,9].apply(lambda x: isinstance(x, float))
review_df_1.loc[mask,'review_title'] = 'No_review_title'

#skin_tone
mask = review_df_1.iloc[:,10].apply(lambda x: isinstance(x, float))
review_df_1.loc[mask,'skin_tone'] = 'No_profile'

#eye_color
mask = review_df_1.iloc[:,11].apply(lambda x: isinstance(x, float))
review_df_1.loc[mask,'eye_color'] = 'No_profile'

#skin_type
mask = review_df_1.iloc[:,12].apply(lambda x: isinstance(x, float))
review_df_1.loc[mask,'skin_type'] = 'No_profile'

#hair_color
mask = review_df_1.iloc[:,13].apply(lambda x: isinstance(x, float))
review_df_1.loc[mask,'hair_color'] = 'No_profile'

2.1.2 Second Partial Review Dataset

review_df_2 = pd.read_csv("dataset/Sephora/reviews_250-500.csv", index_col=0)
column_type_check(review_df_2)

Columns that have multiple data types: 
  review_text
  review_title
  skin_tone
  eye_color
  skin_type
  hair_color

review_df_2.loc[~ review_df_2['author_id'].str.isnumeric(), 'author_id'].unique()

array(['orderGen1254820', 'orderGen1221842', 'orderGen1698648',
       'orderGen53499', 'orderGen51156', 'orderGen333757',
       'orderGen5563740'], dtype=object)

list(author_id_mapping.keys())

['orderGen51156',
 'orderGen2124216',
 'orderGen703225',
 'orderGen5563740',
 'orderGen270100',
 'orderGen1221842',
 'orderGen1254820',
 'orderGen1253445',
 'orderGen1937304',
 'orderGen3046665',
 'orderGen1711826',
 'orderGen309293',
 'orderGen1698648',
 'orderGen39837',
 'orderGen899861']

# author_id
print('author_id that have been replaced with a new id: ', list(author_id_mapping.keys()), '\n')
print('unique string id in this partial dataframe: ', review_df_2.loc[~ review_df_2['author_id'].str.isnumeric(), 'author_id'].unique()
      , '\n')

check_idx = np.isin(review_df_2.loc[~ review_df_2['author_id'].str.isnumeric(), 'author_id'].unique(), 
        list(author_id_mapping.keys()))

print('author_id is in the author id mapping: ', check_idx)

to_be_added = review_df_2.loc[~ review_df_2['author_id'].str.isnumeric(), 'author_id'].unique()[~check_idx]

author_id_mapping_conca = {old: new_id for new_id, old in enumerate(to_be_added, start= 17)}

author_id_mapping = author_id_mapping | author_id_mapping_conca

review_df_2.loc[~ review_df_2['author_id'].str.isnumeric(), 'author_id'] = (
    review_df_2.loc[~ review_df_2['author_id'].str.isnumeric(), 'author_id'].map(author_id_mapping).astype(int)
)

review_df_2['author_id'] = review_df_2['author_id'].astype(int)

author_id that have been replaced with a new id:  ['orderGen51156', 'orderGen2124216', 'orderGen703225', 'orderGen5563740', 'orderGen270100', 'orderGen1221842', 'orderGen1254820', 'orderGen1253445', 'orderGen1937304', 'orderGen3046665', 'orderGen1711826', 'orderGen309293', 'orderGen1698648', 'orderGen39837', 'orderGen899861'] 

unique string id in this partial dataframe:  ['orderGen1254820' 'orderGen1221842' 'orderGen1698648' 'orderGen53499'
 'orderGen51156' 'orderGen333757' 'orderGen5563740'] 

author_id is in the author id mapping:  [ True  True  True False  True False  True]

#review_text
mask = review_df_2.iloc[:,8].apply(lambda x: isinstance(x, float))
review_df_2.loc[mask,'review_text'] = 'No_review'

#review_title
mask = review_df_2.iloc[:,9].apply(lambda x: isinstance(x, float))
review_df_2.loc[mask,'review_title'] = 'No_review_title'

#skin_tone
mask = review_df_2.iloc[:,10].apply(lambda x: isinstance(x, float))
review_df_2.loc[mask,'skin_tone'] = 'No_profile'

#eye_color
mask = review_df_2.iloc[:,11].apply(lambda x: isinstance(x, float))
review_df_2.loc[mask,'eye_color'] = 'No_profile'

#skin_type
mask = review_df_2.iloc[:,12].apply(lambda x: isinstance(x, float))
review_df_2.loc[mask,'skin_type'] = 'No_profile'

#hair_color
mask = review_df_2.iloc[:,13].apply(lambda x: isinstance(x, float))
review_df_2.loc[mask,'hair_color'] = 'No_profile'

2.1.3 Third Partial Review Dataset

review_df_3 = pd.read_csv("dataset/Sephora/reviews_500-750.csv", index_col=0)
column_type_check(review_df_3)

Columns that have multiple data types: 
  review_text
  review_title
  skin_tone
  eye_color
  skin_type
  hair_color

# author_id
print('author_id that have been replaced with a new id: ', list(author_id_mapping.keys()), '\n')
print('unique string id in this partial dataframe: ', review_df_3.loc[~ review_df_3['author_id'].str.isnumeric(), 'author_id'].unique()
      , '\n')

check_idx = np.isin(review_df_3.loc[~ review_df_3['author_id'].str.isnumeric(), 'author_id'].unique(), 
        list(author_id_mapping.keys()))

print('author_id is in the author id mapping: ', check_idx)

to_be_added = review_df_3.loc[~ review_df_3['author_id'].str.isnumeric(), 'author_id'].unique()[~check_idx]

author_id_mapping_conca = {old: new_id for new_id, old in enumerate(to_be_added, start= 19)}

author_id_mapping = author_id_mapping | author_id_mapping_conca

review_df_3.loc[~ review_df_3['author_id'].str.isnumeric(), 'author_id'] = (
    review_df_3.loc[~ review_df_3['author_id'].str.isnumeric(), 'author_id'].map(author_id_mapping).astype(int)
)

review_df_3['author_id'] = review_df_3['author_id'].astype(int)

author_id that have been replaced with a new id:  ['orderGen51156', 'orderGen2124216', 'orderGen703225', 'orderGen5563740', 'orderGen270100', 'orderGen1221842', 'orderGen1254820', 'orderGen1253445', 'orderGen1937304', 'orderGen3046665', 'orderGen1711826', 'orderGen309293', 'orderGen1698648', 'orderGen39837', 'orderGen899861', 'orderGen53499', 'orderGen333757'] 

unique string id in this partial dataframe:  ['orderGen5563740' 'orderGen1474435' 'orderGen1698648'] 

author_id is in the author id mapping:  [ True False  True]

#review_text
mask = review_df_3.iloc[:,8].apply(lambda x: isinstance(x, float))
review_df_3.loc[mask,'review_text'] = 'No_review'

#review_title
mask = review_df_3.iloc[:,9].apply(lambda x: isinstance(x, float))
review_df_3.loc[mask,'review_title'] = 'No_review_title'

#skin_tone
mask = review_df_3.iloc[:,10].apply(lambda x: isinstance(x, float))
review_df_3.loc[mask,'skin_tone'] = 'No_profile'

#eye_color
mask = review_df_3.iloc[:,11].apply(lambda x: isinstance(x, float))
review_df_3.loc[mask,'eye_color'] = 'No_profile'

#skin_type
mask = review_df_3.iloc[:,12].apply(lambda x: isinstance(x, float))
review_df_3.loc[mask,'skin_type'] = 'No_profile'

#hair_color
mask = review_df_3.iloc[:,13].apply(lambda x: isinstance(x, float))
review_df_3.loc[mask,'hair_color'] = 'No_profile'

2.1.4 Fourth Partial Review Dataset

review_df_4 = pd.read_csv("dataset/Sephora/reviews_750-1250.csv", index_col=0)
column_type_check(review_df_4)

C:\Users\fuyuz\AppData\Local\Temp\ipykernel_34616\136717371.py:1: DtypeWarning:

Columns (1) have mixed types. Specify dtype option on import or set low_memory=False.

Columns that have multiple data types: 
  author_id
  review_text
  review_title
  skin_tone
  eye_color
  skin_type
  hair_color

# author_id
mask = review_df_4['author_id'].apply(lambda x: isinstance(x, str))
a = review_df_4.loc[mask, 'author_id']

print('author_id that have been replaced with a new id: ', list(author_id_mapping.keys()), '\n')
print('unique string id in this partial dataframe: ', a[~ a.str.isnumeric()].unique()
      , '\n')

check_idx = np.isin(a[~a.str.isnumeric()].unique(), 
        list(author_id_mapping.keys()))

print('author_id is in the author id mapping: ', check_idx)

to_be_added = a[~a.str.isnumeric()].unique()[~check_idx]

author_id_mapping_conca = {old: new_id for new_id, old in enumerate(to_be_added, start= 20)}

author_id_mapping = author_id_mapping | author_id_mapping_conca


mask_idx = a[a.str.contains(r'order', regex=False)].index # for the str 'order...'

review_df_4.loc[mask_idx, 'author_id'] = review_df_4.loc[mask_idx, 'author_id'].map(author_id_mapping).astype(int)

review_df_4['author_id'] = review_df_4['author_id'].astype(int)

author_id that have been replaced with a new id:  ['orderGen51156', 'orderGen2124216', 'orderGen703225', 'orderGen5563740', 'orderGen270100', 'orderGen1221842', 'orderGen1254820', 'orderGen1253445', 'orderGen1937304', 'orderGen3046665', 'orderGen1711826', 'orderGen309293', 'orderGen1698648', 'orderGen39837', 'orderGen899861', 'orderGen53499', 'orderGen333757', 'orderGen1474435'] 

unique string id in this partial dataframe:  ['orderGen1698648' 'orderGen1566769' 'orderGen1221842' 'orderGen2124216'
 'orderGen3046665' 'orderGen51156' 'orderGen1947347'] 

author_id is in the author id mapping:  [ True False  True  True  True  True False]

#review_text
mask = review_df_4.iloc[:,8].apply(lambda x: isinstance(x, float))
review_df_4.loc[mask,'review_text'] = 'No_review'

#review_title
mask = review_df_4.iloc[:,9].apply(lambda x: isinstance(x, float))
review_df_4.loc[mask,'review_title'] = 'No_review_title'

#skin_tone
mask = review_df_4.iloc[:,10].apply(lambda x: isinstance(x, float))
review_df_4.loc[mask,'skin_tone'] = 'No_profile'

#eye_color
mask = review_df_4.iloc[:,11].apply(lambda x: isinstance(x, float))
review_df_4.loc[mask,'eye_color'] = 'No_profile'

#skin_type
mask = review_df_4.iloc[:,12].apply(lambda x: isinstance(x, float))
review_df_4.loc[mask,'skin_type'] = 'No_profile'

#hair_color
mask = review_df_4.iloc[:,13].apply(lambda x: isinstance(x, float))
review_df_4.loc[mask,'hair_color'] = 'No_profile'

2.1.5 Fifth Partial Review Dataset

review_df_5 = pd.read_csv("dataset/Sephora/reviews_1250-end.csv", index_col=0)
column_type_check(review_df_5)

Columns that have multiple data types: 
  author_id
  review_text
  review_title
  skin_tone
  eye_color
  skin_type
  hair_color

C:\Users\fuyuz\AppData\Local\Temp\ipykernel_34616\1811167476.py:1: DtypeWarning:

Columns (1) have mixed types. Specify dtype option on import or set low_memory=False.

# author_id
mask = review_df_5['author_id'].apply(lambda x: isinstance(x, str))
a = review_df_5.loc[mask, 'author_id']

print('author_id that have been replaced with a new id: ', list(author_id_mapping.keys()), '\n')
print('unique string id in this partial dataframe: ', a[~ a.str.isnumeric()].unique()
      , '\n')

check_idx = np.isin(a[~a.str.isnumeric()].unique(), 
        list(author_id_mapping.keys()))

print('author_id is in the author id mapping: ', check_idx)

mask_idx = a[a.str.contains(r'order', regex=False)].index # for the str 'order...'

review_df_5.loc[mask_idx, 'author_id'] = review_df_5.loc[mask_idx, 'author_id'].map(author_id_mapping).astype(int)

review_df_5['author_id'] = review_df_5['author_id'].astype(int)

author_id that have been replaced with a new id:  ['orderGen51156', 'orderGen2124216', 'orderGen703225', 'orderGen5563740', 'orderGen270100', 'orderGen1221842', 'orderGen1254820', 'orderGen1253445', 'orderGen1937304', 'orderGen3046665', 'orderGen1711826', 'orderGen309293', 'orderGen1698648', 'orderGen39837', 'orderGen899861', 'orderGen53499', 'orderGen333757', 'orderGen1474435', 'orderGen1566769', 'orderGen1947347'] 

unique string id in this partial dataframe:  ['orderGen1947347' 'orderGen1698648' 'orderGen3046665'] 

author_id is in the author id mapping:  [ True  True  True]

#review_text
mask = review_df_5.iloc[:,8].apply(lambda x: isinstance(x, float))
review_df_5.loc[mask,'review_text'] = 'No_review'

#review_title
mask = review_df_5.iloc[:,9].apply(lambda x: isinstance(x, float))
review_df_5.loc[mask,'review_title'] = 'No_review_title'

#skin_tone
mask = review_df_5.iloc[:,10].apply(lambda x: isinstance(x, float))
review_df_5.loc[mask,'skin_tone'] = 'No_profile'

#eye_color
mask = review_df_5.iloc[:,11].apply(lambda x: isinstance(x, float))
review_df_5.loc[mask,'eye_color'] = 'No_profile'

#skin_type
mask = review_df_5.iloc[:,12].apply(lambda x: isinstance(x, float))
review_df_5.loc[mask,'skin_type'] = 'No_profile'

#hair_color
mask = review_df_5.iloc[:,13].apply(lambda x: isinstance(x, float))
review_df_5.loc[mask,'hair_color'] = 'No_profile'

2.1.6 Combine All Datasets

After all the corrections, we combine all our partial review datasets.

Additionally, we assign all the texts a unique ID.

df_lis = [review_df_1, review_df_2, review_df_3, review_df_4, review_df_5]

review_df_all = pd.concat(df_lis, ignore_index=True)

# assign each review a unique ID
review_df_all['review_id'] = pd.factorize(review_df_all['review_text'])[0]

review_df_all.head()

	author_id	rating	is_recommended	helpfulness	total_feedback_count	total_neg_feedback_count	total_pos_feedback_count	submission_time	review_text	review_title	skin_tone	eye_color	skin_type	hair_color	product_id	product_name	brand_name	price_usd	review_id
0	1741593524	5	1.0	1.0	2	0	2	2023-02-01	I use this with the Nudestix “Citrus Clean Bal...	Taught me how to double cleanse!	No_profile	brown	dry	black	P504322	Gentle Hydra-Gel Face Cleanser	NUDESTIX	19.0	0
1	31423088263	1	0.0	NaN	0	0	0	2023-03-21	I bought this lip mask after reading the revie...	Disappointed	No_profile	No_profile	No_profile	No_profile	P420652	Lip Sleeping Mask Intense Hydration with Vitam...	LANEIGE	24.0	1
2	5061282401	5	1.0	NaN	0	0	0	2023-03-21	My review title says it all! I get so excited ...	New Favorite Routine	light	brown	dry	blonde	P420652	Lip Sleeping Mask Intense Hydration with Vitam...	LANEIGE	24.0	2
3	6083038851	5	1.0	NaN	0	0	0	2023-03-20	I’ve always loved this formula for a long time...	Can't go wrong with any of them	No_profile	brown	combination	black	P420652	Lip Sleeping Mask Intense Hydration with Vitam...	LANEIGE	24.0	3
4	47056667835	5	1.0	NaN	0	0	0	2023-03-20	If you have dry cracked lips, this is a must h...	A must have !!!	light	hazel	combination	No_profile	P420652	Lip Sleeping Mask Intense Hydration with Vitam...	LANEIGE	24.0	4

2.2 Review Duplicate

After assigning a unique text ID, we notice that a single product might have multiple repeated reviews. Below are some examples:

mask = review_df_all['review_id'].value_counts()>1
duplicate_idx = review_df_all['review_id'].value_counts()[mask].index
duplicate_idx = duplicate_idx[1:] # filter NA (entried as 'No_review')

duplcate_df = review_df_all[review_df_all['review_id'].isin(duplicate_idx)]
print(duplcate_df.groupby(['product_id', 'review_id'])['review_id'].agg(['count']).sort_values(by = 'count', ascending=False).head())

# non-duplicate dataframe
non_duplicate_idx = review_df_all['review_id'].value_counts()[~mask].index
#non_duplcate_df = review_df_all[review_df_all['review_id'].isin(duplicate_idx)]

                      count
product_id review_id       
P377368    782703        59
P384537    730674        14
P139000    292121        13
P122661    619906        13
P384537    730673         9

What we want is a single product ID corresponds to a unqiue review ID. Let’s remove those duplicates. And our strategy is to keep the first occurence.

mask = review_df_all["review_id"].isin(duplicate_idx)
review_df_all_deduplicate = review_df_all.loc[mask].drop_duplicates(subset=['product_id', 'review_id'], keep='first')

Check again.

print('de-duplicates data frame: ', '\n',
      review_df_all_deduplicate.groupby(['product_id', 'review_id'])['review_id'].agg(['count']).value_counts()
)

de-duplicates data frame:  
 count
1        241827
Name: count, dtype: int64

Combine processed data frame and non-duplicate data frame. Some indixes of the processed data frame are in the non-duplicate data frame.

print('number of processed data frame indixes that are in the non-duplicate data frame: ', np.isin(review_df_all_deduplicate.index,non_duplicate_idx).sum())

number of processed data frame indixes that are in the non-duplicate data frame:  199565

# only need add processed idx that are not in the non-duplicate data frame
idx_to_be_added = review_df_all_deduplicate.index[~ np.isin(review_df_all_deduplicate.index,non_duplicate_idx)]
non_duplicate_idx = np.concatenate([non_duplicate_idx, idx_to_be_added])

review_df_all_deduplicate = review_df_all.loc[non_duplicate_idx]

2.3 NA Value

Check any empty text entries, which we have labeled them as ‘No_review’.

print('number of review entered as \'No_review\': ', (review_df_all_deduplicate['review_text'] == 'No_review').sum())

number of review entered as 'No_review':  1132

review_df_all_deduplicate = review_df_all_deduplicate[review_df_all_deduplicate['review_text'] != 'No_review']

3. EDA

Now we can do the EDA to better understand our data.

3.1 Product Level

product_df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 8494 entries, 0 to 8493
Data columns (total 27 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   product_id          8494 non-null   object 
 1   product_name        8494 non-null   object 
 2   brand_id            8494 non-null   int64  
 3   brand_name          8494 non-null   object 
 4   loves_count         8494 non-null   int64  
 5   rating              8216 non-null   float64
 6   reviews             8216 non-null   float64
 7   size                6863 non-null   object 
 8   variation_type      7050 non-null   object 
 9   variation_value     6896 non-null   object 
 10  variation_desc      1250 non-null   object 
 11  ingredients         7549 non-null   object 
 12  price_usd           8494 non-null   float64
 13  value_price_usd     451 non-null    float64
 14  sale_price_usd      270 non-null    float64
 15  limited_edition     8494 non-null   int64  
 16  new                 8494 non-null   int64  
 17  online_only         8494 non-null   int64  
 18  out_of_stock        8494 non-null   int64  
 19  sephora_exclusive   8494 non-null   int64  
 20  highlights          6287 non-null   object 
 21  primary_category    8494 non-null   object 
 22  secondary_category  8486 non-null   object 
 23  tertiary_category   7504 non-null   object 
 24  child_count         8494 non-null   int64  
 25  child_max_price     2754 non-null   float64
 26  child_min_price     2754 non-null   float64
dtypes: float64(7), int64(8), object(12)
memory usage: 1.7+ MB

First we want to know:

How many brands in this product dataset? How many products in this product dataset?

print('Brand number: ', len(product_df['brand_name'].unique()))
print('Product number: ', len(product_df['product_name'].unique()))

Brand number:  304
Product number:  8415

3.1.1 Rating

Now, let’s explore the average rating of the products. Notice that the value of rating is continuously distributed in [0, 5]. We discretize the variable first.

product_df['rating'].head()

0    3.6364
1    4.1538
2    4.2500
3    4.4762
4    3.2308
Name: rating, dtype: float64

Plot distribution.

bins = [0, 1, 2, 3, 4, 5]
labels = ["(0,1]", "(1,2]", "(2,3]", "(3,4]", "(4,5]"]
product_df["rating_interval"] = pd.cut(product_df["rating"], bins=bins, right=True, labels=labels, include_lowest=False)

ax = sns.barplot(product_df["rating_interval"].value_counts().sort_index())
ax.bar_label(ax.containers[0], fontsize = 10)
ax.set_title('Distribution of Product Rating')
plt.grid(True)
plt.show()

Most of the products are at or above the average level. Let’s also draw a suburst chart to explore the relation between rating and First product category, Second product category , Third product category.

df_plot = product_df[~ product_df['tertiary_category'].isna()] # the chart is not robust to na value

print('Number of rating NA value: ', product_df['tertiary_category'].isna().sum())

fig = px.sunburst(
    df_plot,
    path = ['primary_category', 'secondary_category', 'tertiary_category'],
    color = 'rating',
    maxdepth=-1,
    labels={'rating': 'Product Rating'}
)

fig.update_layout(
    width = 1200,
    height = 900,
    title={
        'text':"Proudct Category and Rating",
        'x': 0.47,
        'y': 0.95,
        'xanchor': 'center'
    },
    margin=dict(t=70)
)

fig.show()
#fig.write_html("my_sunburst_chart.html", include_plotlyjs=True)

Number of rating NA value:  990

Many sub-categories under Fragrance section show deep color, indicating significantly lower ratings compared to the bright orange/yellow. The users are more easily dissatisfied with Sephora’s fragrance-type product offerings.

Skincare (such as “Moisturizers” and “Face Oils”) is dominated by bright yellow and light orange hues, suggesting consistently high customer satisfaction. Makeup, fragrance, and hair show darker orange and reddish tones. This indicates while these categories have a massive product volume, it is more likely to receive lower ratings, whereas skincare serves as the reliable, high-quality backbone of the catalog.

After exploring the relation between product rating and category, how about the product price and product categories?

3.1.2 Price

fig = px.sunburst(
    df_plot,
    path = ['primary_category', 'secondary_category', 'tertiary_category'],
    color = 'price_usd',
    maxdepth=-1,
    labels={'price_usd': 'Product Price'}
)

fig.update_layout(
    width = 1200,
    height = 900,
    title={
        'text':"Proudct Category and Price",
        'x': 0.47,
        'y': 0.95,
        'xanchor': 'center'
    },
    margin=dict(t=70)
)

fig.show()

Notice that there is a very thin but distinct labeled “High Tech Tools” (and nearby “Anti-Aging”) that appears bright yellow, this indicates while the majority of the category is affordable, this specific sub-category represents the premium price ceiling of Sephora’s inventory, standing out sharply against the rest of the dark purple chart.

The massive makeup section (eye, face, lip …) is almost entirely deep purple, indicating a consistently low price point compared to other categories. Makeup unlike skincare or fragrance which have more price variation, has a low barrier to entry cost-wise, it is an accessible category.

By contrast, the overall fragrance section is colored in a lighter tone, indicating that the fragrance products generally sit in a “mid-to-high” price tier.

How about we summarise all the numerical variables in the product dataset?

3.1.3 All Numerical Features

df_plot = product_df[~product_df['reviews'].isna()]

fig = px.parallel_coordinates(
    df_plot.iloc[:,[4, 5, 6, 12]], 
    color="price_usd", 
    color_continuous_scale=px.colors.diverging.Spectral, 
    labels={'price_usd': 'Product Price'},
)

fig.update_layout(
    title={
        'text': "Product Price against Product Numeric Features",
        'x': 0.5,
        'y': 0.99,
        'xanchor': 'center'
    },
    margin=dict(t=80)
)

fig.show()

As the price peaks at the far right, the line dips to near zero on both the reviews and loves_count axes (blue line). There is a clear trade-off between price and engagement. Premium products (e.g. “high tech tools” category) generate almost no community buzz compared to the rest of the categories, likely due to low sales volume.

We see that the red lines (low price products) dominate the top peaks of the loves_count (reaching 1.4M+) and the reviews (reaching 20k+). This shows customers are in favor of cheaper products, the virality is exclusive to the lower price tier.

How about the non-numeric feature?

3.1.4 highlights

Highlights of a product recorded is a list of tags or features that highlight the product’s attributes (e.g. [‘Vegan’, ‘Matte Finish’]).

We extract the top 20 most common keywords across all the highlights of products. Then compare their rating, loves_count, and reviews.

This envolves the creating of uni-gram. We will write a function for creating n-grams.

Note: before generating the unigram, cleaned text is required, so we will also operate necessary text cleaning for the texts in highlights column.

# create n gram
def generate_ngrams(text, n_gram=1):
    token = [token for token in text.lower().split(' ') if token != '' if token not in STOPWORDS]
    ngrams = zip(*[token[i:] for i in range(n_gram)])
    return [' '.join(ngram) for ngram in ngrams]

df_highlights = product_df[~ product_df['highlights'].isna()].copy()

df_highlights['highlights_list'] = df_highlights['highlights'].apply(ast.literal_eval)
df_highlights['highlights_list'] = df_highlights['highlights_list'].apply(lambda x: ' '.join(x))

df_highlights['highlights_list'] = df_highlights['highlights_list'].str.replace(r'(?<!\w)-(?!\w)', ' ', regex=True)
df_highlights['highlights_list'] = df_highlights['highlights_list'].str.replace(r'[^\w\s-]', ' ', regex=True)
df_highlights['highlights_list'] = df_highlights['highlights_list'].str.strip()
df_highlights['highlights_list'] = df_highlights['highlights_list'].str.split()
df_highlights['highlights_list'] = df_highlights['highlights_list'].apply(lambda x: ' '.join(x))

df_highlights['unigram'] = df_highlights['highlights_list'].apply(lambda x: generate_ngrams(x, n_gram=1))

unigrams = defaultdict(int)

for row in df_highlights['unigram']:
    for word in row:
        unigrams[word] += 1

df_unigram = pd.DataFrame(sorted(unigrams.items(), key=lambda x: x[1])[::-1])

Plot the corresponding boxplots.

Note: for reviews and loves_count, we apply \(x \rightarrow log(1+x)\) transformation, otherwise the plots will be too skewed.

rows = list(df_unigram.loc[df_unigram.iloc[:, 1] >= 900, 0]) # Top 20
col_names = ["rating", "reviews", "loves_count"]

for col in col_names:
    plot_data = []
    for word in rows:
        sub = df_highlights[df_highlights["highlights_list"]
                            .str.contains(word, case=False, regex=False)][col]
        vals = np.log1p(sub) if col != "rating" else sub
        plot_data.append(pd.DataFrame({
            "word": word,
            col: vals
        }))
    plot_df = pd.concat(plot_data, ignore_index=True)

    plt.figure(figsize=(12, 6))
    sns.boxplot(data=plot_df, x="word", y=col, color="skyblue")
    plt.xticks(rotation=45, ha="right")
    plt.xlabel("keyword")
    plt.ylabel(col if col == "rating" else f"log (1 + {col})")
    plt.title(f"{col} by Top-20 highlights keyword")
    plt.tight_layout()
    plt.show()

Based on the 3 plots above, there is no much variation across the keywords, highlight keyword is not the sole factor impacting the numercial metrics here.

3.2 Review Level

First, how many customers posted at least 1 review?

print('Unique customers: ', len(review_df_all_deduplicate['author_id'].unique()))

Unique customers:  432578

Then, what does the time interval between posted reviews look like?

3.2.1 Time Range

review_df_all_deduplicate['submission_year'] = pd.to_datetime(review_df_all_deduplicate['submission_time']).dt.year
review_df_all_deduplicate['submission_month'] = pd.to_datetime(review_df_all_deduplicate['submission_time']).dt.month

plot_data = review_df_all_deduplicate.loc[:, ['submission_year', 'submission_month']]
plot_data = plot_data.groupby(["submission_year", "submission_month"]).size().reset_index(name="count")

fig = px.sunburst(
    plot_data,
    path=["submission_year", "submission_month"],
    values="count",
    color="count",
    color_continuous_scale="YlOrRd",
    labels={"count": "Reviews Count by Month"},
)

fig.update_layout(
    width = 1200,
    height = 900,
    title={
        'text':"Reviews Count",
        'x': 0.47,
        'y': 0.95,
        'xanchor': 'center'
    },
    margin=dict(t=70)
)

fig.show()

Note: ignore the color presentation (difference) in the middle pie chart.

The reviews span from 2008 to 2023. Among these years, 2020, 2021, and 2022 account for a large proportion of the reviews. In addition, reviews peak in January, April, May, and August during these years.

Since we are dealing with unstructured text data, N-gram analysis is essential to capture recurring phrase patterns and understand the linguistic structure of the reviews.

3.2.2 N-Gram Analysis

Create dataframe containing n-grams needed, we will be exploring unigram, bigram, and trigram.

df_gram = review_df_all_deduplicate.copy()
df_gram = df_gram[['author_id','review_id','review_text']]

# unigram
df_gram['unigram'] = review_df_all_deduplicate['review_text'].apply(lambda x: generate_ngrams(x, n_gram=1))
# bigram
df_gram['bigram'] = review_df_all_deduplicate['review_text'].apply(lambda x: generate_ngrams(x, n_gram=2))
# trigram
df_gram['trigram'] = review_df_all_deduplicate['review_text'].apply(lambda x: generate_ngrams(x, n_gram=3))

Unigram

Most common unigrams are mostly stop words and uncleand words (with punctuations), which do not contain much information, this indicates we need further cleaning.

unigrams = defaultdict(int)

for row in df_gram['unigram']:
    for word in row:
        unigrams[word] += 1

df_unigram = pd.DataFrame(sorted(unigrams.items(), key=lambda x: x[1])[::-1])

fig, ax = plt.subplots(figsize=(18, 50), dpi=100)
N = 100
sns.barplot(y = df_unigram[0][:N], x=df_unigram[1][:N])

ax.set_xlabel('')
ax.set_ylabel('')
ax.set_title('Top 100 most common unigrams in customer reviews')
plt.show()

Bigram

Bigrams reveal more information than the unigrams, the most common bigrams are about skin, deatiled information about customer skin is shown, for instance, we see ‘dry skin’ and ‘sensitive skin’ are mentioned very frequently. And we can infer that people who bought skincare products are likely to leave a comment.

But most bigrams still contain stop words and uncleand words (with punctuations), we need further cleaning.

bigrams = defaultdict(int)

for row in df_gram['bigram']:
    for word in row:
        bigrams[word] += 1

df_bigram = pd.DataFrame(sorted(bigrams.items(), key=lambda x: x[1])[::-1])

fig, ax = plt.subplots(figsize=(18, 50), dpi=100)

sns.barplot(y = df_bigram[0][:N], x=df_bigram[1][:N])

ax.set_xlabel('')
ax.set_ylabel('')
ax.set_title('Top 100 most common bigrams in customer reviews')
plt.show()

Trigram

Most common trigrams are also about skin, from the most common trigram we see customers care long-lasting effects of the prodcuts.

The same stop words and uncleand words issues exist in this plot, we need further cleaning.

trigrams = defaultdict(int)

for row in df_gram['trigram']:
    for word in row:
        trigrams[word] += 1

df_trigram = pd.DataFrame(sorted(trigrams.items(), key=lambda x: x[1])[::-1])

fig, ax = plt.subplots(figsize=(18, 50), dpi=100)

sns.barplot(y = df_trigram[0][:N], x=df_trigram[1][:N])

ax.set_xlabel('')
ax.set_ylabel('')
ax.set_title('Top 100 most common trigrams in customer reviews')
plt.show()

This is enough for EDA, now as the above N-Gram analysis indicated, we proceed to clean texts before word embedding.

4. Text Cleaning

4.1 Language Detection

This dataset contains multiple language product reviews, examples are Spanish, French, Chinese, and etc., so our first thing to do is to detect English-only reviews.

We will use two language detection tools: langdetect and langid, the first one will do the detection work, and the second one will validate the detection results. Detailed scripts are omitted for brevity, interesting reader can find it here.

Here we directly import the processed dataset.

review_df_all_deduplicate = pd.read_csv("dataset/review_df_all_deduplicate_english.csv", index_col = 0)

4.2 Lowercase Letter

review_df_all_deduplicate['review_text'] = review_df_all_deduplicate['review_text'].str.lower()

4.3 Standardize Unicode Character

After detecting the non-English reviews, there are other unicode variation of standard punctuation marks,e.g. ‘！’ \(\rightarrow\) ‘!’, we clean them by mapping them to standard English punctuation characters.

We first use unicodedata library to systematically normalize the variation, then utilize LLM to detect remaining variation and list the corresponding character mapping, the detail of using LLM is ommited.

def clean_unicode(text):
    text = unicodedata.normalize('NFKC', text)

    return text

def detect_anomalies(text):
    anomalies = []
    for char in text:
        cp = ord(char)
        # unicode point for ASCII (0-127) 
        if cp < 128:
            continue
            
        category = unicodedata.category(char)
        # category of unicode: Punctuation, other (PO)
        if category in ['Po']:
            anomalies.append(f"{char}") # ({hex(cp)})
            
    return ",".join(anomalies) if anomalies else None

review_df_all_deduplicate['review_text_cleaned'] = review_df_all_deduplicate['review_text'].map(clean_unicode)

a = review_df_all_deduplicate['review_text_cleaned'].apply(detect_anomalies)
# export for LLM to detect
a[a.notna()].to_csv('dataset/special_punctuation.csv', encoding="utf-8-sig", index=True)

def clean_unicode(text):
    punctuation_map = {
        '。': '.',        
        '、': ',',   
        '·': '.',         
        '・': '.',       
        '¡': '!',      
        '¿': '?',        
        '،': ',',  
        }
    
    unmappable_symbols = ['•', '⁃']

    for old, new in punctuation_map.items():
        text = re.sub(old, new, text)

    for char in unmappable_symbols:
        text = re.sub(char, ' ', text)

    return text

review_df_all_deduplicate['review_text_cleaned'] = review_df_all_deduplicate['review_text_cleaned'].map(clean_unicode)

4.4 Non-ASCII Character

This is similar to unicode cleaning, here we will remove those unmeaningful non-ASCII characters.

Non-ASCII characters are such as non-English characters ‘å’, emojis ‘🙂’, and the characters that cannot be displayed ‘028’. However, we will only remain emojis, since the later embedding model can identify and convert emojis to correct numerical vectors.

# English letters + number + ASCII punctuation
ascii_allowed = r"[A-Za-z0-9\s\.,!?;:'\"()\-\[\]{}]"

emoji_pattern = (
    r"[\U0001F1E0-\U0001F1FF"  # Flags
    r"\U0001F300-\U0001F5FF"   # Symbols & pictographs
    r"\U0001F600-\U0001F64F"   # Emoticons
    r"\U0001F680-\U0001F6FF"   # Transport & map symbols
    r"\U0001F700-\U0001F77F"   # Alchemical symbols
    r"\U0001F780-\U0001F7FF"   # Geometric symbols
    r"\U0001F800-\U0001F8FF"   # Supplemental arrows
    r"\U0001F900-\U0001F9FF"   # Supplemental symbols & pictographs
    r"\U0001FA00-\U0001FA6F"   # Chess, symbols
    r"\U0001FA70-\U0001FAFF"   # Emoji components
    r"\U00002702-\U000027B0"   # Dingbats
    r"\U000024C2-\U0001F251"   # Enclosed characters
 r"]"
)

allowed = f"(?:{ascii_allowed}|{emoji_pattern})"

def clean_text_keep_emoji(text):
    return "".join(char if re.match(allowed, char) else "" for char in text)

review_df_all_deduplicate['review_text_cleaned'] = review_df_all_deduplicate['review_text_cleaned'].apply(lambda x: clean_text_keep_emoji(x))

punc_lis_remain = []

for text in review_df_all_deduplicate['review_text_cleaned']:
    punc_lis_remain.extend(re.findall(r"[^\w\s]", text))

punc_lis_remain = sorted(set(punc_lis_remain), key=lambda ch: ord(ch))
punc_lis_remain = [p for p in punc_lis_remain if not p.isascii()]

print('Non-ASCII characters remained: ', punc_lis_remain)

Non-ASCII characters remained:  ['─', '╥', '■', '▪', '▫', '▰', '▶', '◇', '◡', '◼', '☀', '☁', '★', '☆', '☝', '☬', '☹', '☺', '♀', '♂', '♡', '♥', '♻', '♾', '⚗', '⚜', '⚠', '⚡', '⚪', '⛅', '⛑', '✅', '✈', '✋', '✌', '✓', '✔', '✖', '✦', '✨', '✳', '❄', '❌', '❕', '❗', '❣', '❤', '➕', '➖', '➡', '➾', '➿', '⠀', '⬇', '⭐', '〈', '《', '》', '〜', '〰', '\uf04c', '\uf0a7', '︅', '︎', '️', '\ufeff', '', '�', '🅶', '🅷', '🅼', '🅾', '🆂', '🆃', '🆄', '🆈', '🆓', '🆙', '🌊', '🌙', '🌞', '🌟', '🌫', '🌱', '🌸', '🌹', '🌺', '🌿', '🍂', '🍃', '🍅', '🍉', '🍊', '🍋', '🍌', '🍒', '🍓', '🍞', '🍬', '🍯', '🎀', '🎄', '🎉', '🏆', '🏻', '🏼', '🏽', '🏾', '🐐', '🐘', '🐝', '🐣', '👀', '👁', '👄', '👋', '👌', '👍', '👎', '👏', '👣', '👩', '👶', '💀', '💄', '💅', '💋', '💓', '💔', '💕', '💖', '💗', '💘', '💙', '💚', '💛', '💜', '💞', '💡', '💥', '💦', '💨', '💫', '💯', '📌', '📍', '📦', '🔥', '🔵', '🖤', '😀', '😁', '😂', '😃', '😄', '😅', '😆', '😉', '😊', '😌', '😍', '😏', '😒', '😓', '😔', '😕', '😖', '😘', '😚', '😝', '😞', '😠', '😢', '😩', '😪', '😫', '😬', '😭', '😮', '😲', '😳', '😶', '😻', '🙁', '🙂', '🙃', '🙄', '🙌', '🙏', '🛍', '🛑', '🛒', '🤌', '🤍', '🤎', '🤐', '🤓', '🤔', '🤗', '🤞', '🤡', '🤢', '🤣', '🤤', '🤦', '🤨', '🤩', '🤪', '🤭', '🤮', '🤯', '🤷', '🥑', '🥰', '🥲', '🥳', '🥴', '🥶', '🥹', '🥺', '🦝', '🦲', '🧖', '🧡', '🧴', '🧼', '🧿', '🩸', '🪄', '🫠', '🫣', '🫰', '🫶']

non_emoji_list = ['─', '╥', '■', '▪', '▫', '▰', '▶', '◇', '◡', '◼', 
                  '〈', '《', '》', '〜', '〰', 
                  '✦', '✳', '⠀', '︅', '︎', '️', '\ufeff', '', '\uf04c', '\uf0a7']

pat = "[" + re.escape("".join(non_emoji_list)) + "]"
review_df_all_deduplicate["review_text_cleaned"] = review_df_all_deduplicate["review_text_cleaned"].str.replace(pat, "", regex=True)

punc_lis_remain = []
for text in review_df_all_deduplicate['review_text_cleaned']:
    punc_lis_remain.extend(re.findall(r"[^\w\s]", text))

punc_lis_remain = sorted(set(punc_lis_remain), key=lambda ch: ord(ch))
punc_lis_remain = [p for p in punc_lis_remain if not p.isascii()]

print('Non-ASCII characters remained after cleaning: ', punc_lis_remain)

Non-ASCII characters remained after cleaning:  ['☀', '☁', '★', '☆', '☝', '☬', '☹', '☺', '♀', '♂', '♡', '♥', '♻', '♾', '⚗', '⚜', '⚠', '⚡', '⚪', '⛅', '⛑', '✅', '✈', '✋', '✌', '✓', '✔', '✖', '✨', '❄', '❌', '❕', '❗', '❣', '❤', '➕', '➖', '➡', '➾', '➿', '⬇', '⭐', '�', '🅶', '🅷', '🅼', '🅾', '🆂', '🆃', '🆄', '🆈', '🆓', '🆙', '🌊', '🌙', '🌞', '🌟', '🌫', '🌱', '🌸', '🌹', '🌺', '🌿', '🍂', '🍃', '🍅', '🍉', '🍊', '🍋', '🍌', '🍒', '🍓', '🍞', '🍬', '🍯', '🎀', '🎄', '🎉', '🏆', '🏻', '🏼', '🏽', '🏾', '🐐', '🐘', '🐝', '🐣', '👀', '👁', '👄', '👋', '👌', '👍', '👎', '👏', '👣', '👩', '👶', '💀', '💄', '💅', '💋', '💓', '💔', '💕', '💖', '💗', '💘', '💙', '💚', '💛', '💜', '💞', '💡', '💥', '💦', '💨', '💫', '💯', '📌', '📍', '📦', '🔥', '🔵', '🖤', '😀', '😁', '😂', '😃', '😄', '😅', '😆', '😉', '😊', '😌', '😍', '😏', '😒', '😓', '😔', '😕', '😖', '😘', '😚', '😝', '😞', '😠', '😢', '😩', '😪', '😫', '😬', '😭', '😮', '😲', '😳', '😶', '😻', '🙁', '🙂', '🙃', '🙄', '🙌', '🙏', '🛍', '🛑', '🛒', '🤌', '🤍', '🤎', '🤐', '🤓', '🤔', '🤗', '🤞', '🤡', '🤢', '🤣', '🤤', '🤦', '🤨', '🤩', '🤪', '🤭', '🤮', '🤯', '🤷', '🥑', '🥰', '🥲', '🥳', '🥴', '🥶', '🥹', '🥺', '🦝', '🦲', '🧖', '🧡', '🧴', '🧼', '🧿', '🩸', '🪄', '🫠', '🫣', '🫰', '🫶']

4.5 Abbreviation

Two types of abbreviation:

1. Contraction e.g.: i’ve and don’t, …

2. Abbreviation:

   • slang/Internet slang, e.g.: idk, lol, omg, …

   • specific terms, e.g: aha, bha, …

For type 1:

We detect 0.

# previous n gram generator function, we slightly fix condition, don't let stopwords removed
def generate_ngrams_1(text, n_gram=1):
    token = [token for token in text.lower().split(' ') if token != '']
    ngrams = zip(*[token[i:] for i in range(n_gram)])
    return [' '.join(ngram) for ngram in ngrams]

df_gram_abbreviation = pd.DataFrame(review_df_all_deduplicate['review_text_cleaned'])
df_gram_abbreviation['unigram'] = df_gram_abbreviation['review_text_cleaned'].apply(lambda x: generate_ngrams_1(x, n_gram=1))

all_tokens = (tok for lst in df_gram_abbreviation["unigram"] for tok in lst)

pattern = re.compile(r"[']")
abbr_tokens = [t for t in all_tokens if pattern.search(t)]

print('Number of tokens that contain \': ', len(abbr_tokens))

Number of tokens that contain ':  0

For type 2:

df_gram_abbreviation = pd.DataFrame(review_df_all_deduplicate['review_text_cleaned'])
df_gram_abbreviation['unigram'] = df_gram_abbreviation['review_text_cleaned'].apply(lambda x: generate_ngrams(x, n_gram=1))

all_tokens = (tok for lst in df_gram_abbreviation["unigram"] for tok in lst)
short_tokens = [t for t in all_tokens if len(t) < 4]
short_freq = Counter(short_tokens)

print('Top 300 frequent abbreviation: ')
# top 300
short_freq.most_common(300) 

Top 300 frequent abbreviation: 

[('use', 272304),
 ('ive', 171096),
 ('dry', 159353),
 ('im', 154930),
 ('one', 127205),
 ('it.', 124566),
 ('eye', 78190),
 ('see', 77121),
 ('try', 74632),
 ('day', 67153),
 ('got', 64184),
 ('oil', 59045),
 ('now', 57219),
 ('bit', 56869),
 ('put', 49486),
 ('way', 48134),
 ('say', 44528),
 ('me.', 42047),
 ('-', 41424),
 ('it!', 41128),
 ('lot', 40673),
 ('new', 40599),
 ('lip', 40348),
 ('go', 40282),
 ('buy', 39846),
 ('it,', 34403),
 ('two', 32847),
 ('2', 31272),
 ('3', 25486),
 ('far', 24733),
 ('bad', 22231),
 ('ill', 21496),
 ('id', 19492),
 ('job', 18429),
 ('gel', 18352),
 ('big', 17864),
 ('top', 17829),
 ('.', 17445),
 ('saw', 17114),
 ('may', 16945),
 ('red', 16641),
 ('spf', 16066),
 ('c', 15669),
 ('5', 15594),
 ('fan', 15298),
 ('yet', 15091),
 ('4', 14276),
 ('up.', 13624),
 ('let', 13066),
 ('end', 12823),
 ('add', 12805),
 ('due', 12034),
 ('me,', 11942),
 ('on.', 11419),
 ('!', 11052),
 ('rid', 10846),
 ('sun', 10810),
 ('jar', 10609),
 ('(i', 10440),
 (',', 9992),
 ('run', 9788),
 ('old', 9731),
 ('10', 9385),
 ('rub', 8722),
 ('bed', 8609),
 ('1', 8498),
 ('mix', 7952),
 (':)', 7738),
 ('100', 7611),
 ('ago', 7481),
 ('set', 6567),
 ('is.', 6520),
 ('box', 6383),
 ('on,', 6206),
 ('do.', 6172),
 ('me!', 6158),
 ('in.', 5749),
 ('non', 5487),
 ('6', 5468),
 ('up,', 5288),
 ('30', 5086),
 ('tan', 5078),
 ('dr.', 4841),
 ('de', 4401),
 ('bc', 4292),
 ('ok', 4279),
 ('2-3', 4275),
 ('la', 4225),
 ('to.', 4189),
 ('pay', 4146),
 ('per', 4145),
 ('hot', 4060),
 ('wow', 4003),
 ('sit', 3997),
 ('20', 3982),
 ('so,', 3861),
 ('def', 3856),
 ('7', 3826),
 ('ran', 3777),
 ('oh', 3601),
 ('tea', 3588),
 ('mom', 3573),
 ('t', 3548),
 ('yes', 3426),
 (':(', 3407),
 ('15', 3395),
 ('wet', 3222),
 ('is,', 3218),
 ('ton', 3134),
 ('pad', 3058),
 ('pat', 2932),
 ('age', 2922),
 ('kit', 2788),
 ('8', 2671),
 ('up!', 2638),
 ('in,', 2570),
 ('tag', 2568),
 ('50', 2516),
 ('omg', 2448),
 ('(', 2443),
 ('sad', 2431),
 ('cut', 2425),
 ('!!', 2424),
 ('aha', 2407),
 ('mer', 2406),
 ('--', 2378),
 ('fun', 2378),
 ('ole', 2373),
 ('dab', 2358),
 ('bag', 2350),
 ('pm', 2343),
 ('40', 2328),
 ('12', 2311),
 ('u', 2288),
 ('lol', 2281),
 ('fix', 2269),
 ('low', 2267),
 ('bar', 2263),
 ('pea', 2164),
 ('go.', 2161),
 ('ptr', 2055),
 ('3-4', 2043),
 ('us', 2034),
 ('pop', 2019),
 ('ok.', 2018),
 ('air', 1980),
 ('1-2', 1960),
 ('be.', 1936),
 ('hit', 1925),
 ('...', 1872),
 ('bha', 1846),
 ('spa', 1838),
 ('do,', 1836),
 ('cc', 1804),
 ('so.', 1780),
 ('fit', 1780),
 ('soo', 1775),
 ('dr', 1771),
 ('w', 1755),
 ('fab', 1731),
 ('lid', 1715),
 ('aid', 1698),
 ('oz', 1688),
 ('ask', 1688),
 ('odd', 1646),
 ('tip', 1644),
 ('tub', 1605),
 ('tad', 1584),
 ('bb', 1567),
 ('(im', 1565),
 ('3rd', 1534),
 ('cap', 1533),
 ('con', 1493),
 ('is!', 1491),
 ('jet', 1489),
 ('sat', 1470),
 ('ha', 1452),
 ('(it', 1451),
 ('2nd', 1444),
 ('on!', 1439),
 ('2x', 1434),
 ('30s', 1428),
 ('key', 1427),
 ('mid', 1424),
 ('ok,', 1392),
 ('boy', 1387),
 ('of.', 1372),
 ('soy', 1362),
 ('zit', 1353),
 ('ren', 1352),
 ('(or', 1313),
 ('god', 1303),
 ('45', 1294),
 ('via', 1289),
 ('win', 1287),
 ('lag', 1286),
 ('(my', 1271),
 ('!!!', 1267),
 ('..', 1237),
 ('(as', 1197),
 ('vit', 1195),
 (':', 1190),
 ('24', 1190),
 ('to,', 1176),
 ('tin', 1173),
 ('ur', 1167),
 ('dot', 1158),
 ('(a', 1154),
 ('ten', 1131),
 ('no.', 1129),
 ('60', 1115),
 ('20s', 1103),
 ('pot', 1093),
 ('bye', 1089),
 ('(in', 1082),
 ('idk', 1067),
 ('25', 1065),
 ('n', 1059),
 ('1st', 1058),
 ('fav', 1047),
 ('e', 1045),
 ('min', 1029),
 ('14', 997),
 ('(if', 985),
 ('tho', 983),
 ('(no', 982),
 ('six', 978),
 ('hg', 974),
 ('row', 965),
 ('jaw', 961),
 ('dew', 952),
 ('❤', 949),
 ('c.', 937),
 ('c,', 917),
 ('etc', 911),
 ('die', 894),
 ('bay', 874),
 ('vs', 867),
 ('vib', 864),
 ('do!', 855),
 (')', 854),
 ('4-5', 846),
 ('am.', 845),
 ('duo', 831),
 ('to!', 819),
 ('eat', 803),
 ('2.', 797),
 ('90', 794),
 ('dip', 791),
 ('pro', 781),
 ('man', 780),
 ('re', 774),
 ('9', 752),
 ('jlo', 748),
 ('1.', 745),
 ('lil', 743),
 ('35', 735),
 ('ups', 720),
 ('34', 701),
 ('pm.', 697),
 ('vox', 695),
 ('80', 692),
 ('go!', 692),
 ('son', 691),
 ('st.', 682),
 ('s', 680),
 ('4th', 677),
 ('ph', 676),
 ('tlc', 675),
 ('13', 664),
 ('tap', 644),
 ('yo', 637),
 ('uv', 634),
 ('itd', 629),
 ('cuz', 628),
 ('(at', 623),
 ('28', 618),
 (';)', 617),
 ('40s', 611),
 ('mad', 609),
 ('raw', 609),
 ('it?', 609),
 ('18', 604),
 ('o', 602),
 ('3x', 601),
 ('it)', 593),
 ('23', 578),
 ('me)', 576),
 ('(so', 576),
 ('guy', 563),
 ('no,', 562),
 ('pre', 557),
 ('ceo', 557),
 ('0', 556),
 ('gym', 553),
 ('go,', 550),
 ('oat', 550),
 ('55', 549),
 ('ice', 517),
 ('met', 509)]

We standardized abbreviations by examining the top 300 most frequent tokens with a length of less than 4 characters. While these thresholds (Top-300, Length < 4) are heuristic, this approach effectively captures the majority of high-frequency cases without requiring exhaustive manual review.

Create correction mapping.

correct_map = {
    r"\bive\b": "i have",
    r"\bim\b": "i am",
    r"\bill\b": "i will",
    r"\bspf\b": "sun protection factor",
    r"\bbc\b": "because",
    r"\bdr\.\b": "doctor",
    r"\bdr\b": "doctor",
    r"\bdef\b": "definitely",
    r"\bptr\b": "peter thomas roth",
    r"\bbha\b": "beta hydroxy acid",
    r"\baha\b": "alpha hydroxy acid",
    r"\bbb\b": "beauty balm",
    r"\bha\b": "hyaluronic acid",
    r"\bomg\b": "oh my god",
    r"\blol\b": "laugh out loud",
    r"\bidk\b": "i do not know",
    r"\bur\b": "your",
    r"\btho\b": "though",
    r"\bfav\b": "favorite",
    r"\bceo\b": "sunday riley",
    r"\bhg\b": "holy grail",
    r"\betc\b": "et cetera",
    r"\bvib\b": "very important beauty insider",
    r"\bmin\b": "minutes",
    r"\bph\b": "potential of hydrogen",
    r"\bmeh\b": "eh",
    r"\bfab\b": "first aid beauty",
    r"\buv\b": "ultraviolet",
    r"\bitd\b": "it would",
    r"\bcuz\b": "because",
    r"\bvit\b": "vitamin",
    r"\b30s\b": "30 seconds",
    r"\b20s\b": "20 seconds",
    r"\b40s\b": "40 seconds",
    r"\b1st\b": "first",
    r"\b2nd\b": "second",
    r"\b3rd\b": "third",
    r"\b4th\b": "fourth",
    r"\b3x\b": "3 times",
    r"\blil\b": "little",
    
    # latter added
    r"\bwouldnt\b": "would not",
    r"\bhavec\b": "have",
    r"\bbomb.com\b": "great",
    r"\byrs\b": "years",
    r"\bfyi\b": "for your information",
    r"\brn\b": "right now",
    r"\bsoo\b": "so",
    r"\bgonna\b": "going to",
    r"\byoull\b": "you will",

    r"(?<=\s)\(im(?=\s)": "i am",
    r"(?<=\s)\(if(?=\s)": "if",
    r"(?<=\s)\(no(?=\s)": "no",
    r"(?:(?<=\s)|^)\bit\)(?=\s|$)": "it",
    r"(?:(?<=\s)|^)\bme\)(?=\s|$)": "me",
    r"(?<=\s)\(so(?=\s)": "so",
}

review_df_all_deduplicate['review_text_cleaned'] = (
    review_df_all_deduplicate['review_text_cleaned']
    .replace(correct_map,regex=True)
)

Specifically fix ‘ii’ abbrevation.

‘ii’ could be typo of ‘i’, and also could be the product name ‘sk-ii’.

COMMON_VERBS_AFTER_I = {
    "love","like","think","use","used","feel","felt","have","had","am",
    "would","will","wish","want","tried","try","bought","buy",
    "recommend","ordered","order","see","saw","noticed","notice",
    "find","found","hate","dislike","prefer","need","needed",

    # latter added
    "mean", "only", "also", "heard", 
}

def fix_ii_typos(text):

    tokens = text.split()
    new_tokens = []

    for idx, tok in enumerate(tokens):
        lower_tok = tok.lower()

        if lower_tok == "ii":
            prev_tok = tokens[idx-1].lower() if idx > 0 else ""
            next_tok = tokens[idx+1].lower() if idx+1 < len(tokens) else ""

            if prev_tok.isdigit():
                new_tokens.append(tok)

            elif next_tok in COMMON_VERBS_AFTER_I:

                new_tokens.append("i")
            else:

                new_tokens.append(tok)
        else:
            new_tokens.append(tok)

    return " ".join(new_tokens)

mask = review_df_all_deduplicate['review_text_cleaned'].str.contains(r'\bii\b', regex=True)
review_df_all_deduplicate.loc[mask, 'review_text_cleaned'] = review_df_all_deduplicate.loc[mask, 'review_text_cleaned'].apply(fix_ii_typos)

4.6 URL

Clean urls in reviews.

def clean_URL(text):
    url = re.compile(r'https?://\S+|www\.\S+|\bhttps?\b') # enhanced, remove residual 'http'
    
    return url.sub('',text).strip()

review_df_all_deduplicate['review_text_cleaned'] = review_df_all_deduplicate['review_text_cleaned'].map(clean_URL)
review_df_all_deduplicate['review_text_cleaned'] = review_df_all_deduplicate['review_text_cleaned'].str.replace(r"\b[\w.-]+\.com\b", "", regex=True)

4.7 Letter Repeat

Clean the words that contain multiple repeated letters. (e.g., loooooooove, soooooooo …)

1. reduce any character sequence exceeding a length of 3 down to 2 characters to restore standard spelling.

2. Find and analyze the top 500 most frequent words containing repeated characters to identify and handle remaining irregularities.

mask = review_df_all_deduplicate['review_text_cleaned'].str.contains(r"([A-Za-z])\1{2,}", regex=True)


review_df_all_deduplicate.loc[mask, 'review_text_cleaned'] = (
    review_df_all_deduplicate.loc[mask, 'review_text_cleaned']
    .str.replace(r"([A-Za-z])\1{2,}", r"\1\1", regex=True)
)

C:\Users\fuyuz\AppData\Local\Temp\ipykernel_34616\277250214.py:1: UserWarning:

This pattern is interpreted as a regular expression, and has match groups. To actually get the groups, use str.extract.

df_gram_repeat = review_df_all_deduplicate.loc[mask, 'review_text_cleaned'].copy()
df_gram_repeat = pd.DataFrame(df_gram_repeat)
df_gram_repeat['unigram'] = df_gram_repeat['review_text_cleaned'].apply(lambda x: generate_ngrams(x, n_gram=1))

all_tokens = (tok for lst in df_gram_repeat["unigram"] for tok in lst)

short_tokens = [t for t in all_tokens if re.search(r"([A-Za-z])\1", t)] 
short_freq = Counter(short_tokens)

print('Top 500 frequent words with repeated letter: ')

# top 500
short_freq.most_common(500)

Top 500 frequent words with repeated letter: 

[('soo', 7964),
 ('really', 5194),
 ('will', 4676),
 ('feel', 3779),
 ('little', 3375),
 ('feels', 3206),
 ('good', 2891),
 ('smells', 2337),
 ('recommend', 1877),
 ('look', 1801),
 ('smell', 1734),
 ('feeling', 1667),
 ('see', 1619),
 ('well', 1548),
 ('smooth', 1532),
 ('still', 1418),
 ('need', 1384),
 ('stuff', 1285),
 ('looks', 1284),
 ('looking', 1239),
 ('actually', 1159),
 ('better', 1113),
 ('week', 1050),
 ('lovee', 1028),
 ('apply', 1025),
 ('full', 1019),
 ('pretty', 1004),
 ('difference', 996),
 ('bottle', 951),
 ('free', 930),
 ('loove', 914),
 ('weeks', 813),
 ('literally', 812),
 ('less', 796),
 ('happy', 768),
 ('keep', 755),
 ('usually', 743),
 ('looked', 736),
 ('getting', 697),
 ('small', 690),
 ('redness', 662),
 ('sunscreen', 644),
 ('especially', 601),
 ('all.', 599),
 ('good.', 586),
 ('different', 559),
 ('well.', 528),
 ('tell', 519),
 ('applying', 495),
 ('good!', 481),
 ('took', 461),
 ('finally', 446),
 ('totally', 420),
 ('applied', 402),
 ('loong', 398),
 ('smooth.', 394),
 ('overall', 388),
 ('wayy', 359),
 ('better.', 355),
 ('immediately', 352),
 ('keeps', 342),
 ('putting', 341),
 ('seen', 337),
 ('too.', 320),
 ('off.', 319),
 ('seems', 316),
 ('irritate', 316),
 ('recommended', 298),
 ('summer', 298),
 ('add', 292),
 ('waay', 277),
 ('needed', 267),
 ('recommend!', 260),
 ('all,', 257),
 ('cheeks', 251),
 ('three', 250),
 ('smooth,', 243),
 ('soon', 243),
 ('deep', 242),
 ('smaller', 241),
 ('application', 240),
 ('good,', 238),
 ('well,', 233),
 ('gotten', 231),
 ('smell.', 228),
 ('recommend.', 225),
 ('off,', 224),
 ('stuff.', 218),
 ('added', 218),
 ('normally', 216),
 ('difference.', 216),
 ('obsessed', 216),
 ('issues', 214),
 ('bottle.', 212),
 ('peel', 209),
 ('sleeping', 209),
 ('smoother', 207),
 ('cooling', 205),
 ('suuper', 197),
 ('stopped', 197),
 ('gross', 196),
 ('effective', 193),
 ('matte', 189),
 ('irritated', 188),
 ('smell,', 188),
 ('seeing', 187),
 ('too!', 184),
 ('irritation', 183),
 ('impressed', 182),
 ('seem', 178),
 ('personally', 174),
 ('sleep', 171),
 ('effect', 166),
 ('follow', 163),
 ('stuff!', 161),
 ('green', 161),
 ('cool', 160),
 ('week.', 159),
 ('disappointed', 158),
 ('matter', 153),
 ('glass', 153),
 ('waterproof', 151),
 ('essence', 151),
 ('dryness', 150),
 ('fell', 149),
 ('soothing', 147),
 ('bigger', 146),
 ('clogged', 146),
 ('supposed', 143),
 ('chapped', 142),
 ('sunscreen.', 141),
 ('seemed', 140),
 ('longg', 138),
 ('supple', 135),
 ('well!', 135),
 ('overall,', 135),
 ('smooth!', 134),
 ('weeks.', 133),
 ('basically', 133),
 ('guess', 132),
 ('issue', 129),
 ('all!', 128),
 ('peeling', 127),
 ('needs', 125),
 ('typically', 124),
 ('application.', 123),
 ('adding', 121),
 ('nighttime', 121),
 ('worried', 121),
 ('cotton', 121),
 ('keeping', 119),
 ('horrible', 118),
 ('redness,', 117),
 ('rubbing', 117),
 ('immediate', 116),
 ('suffer', 116),
 ('unless', 116),
 ('week,', 115),
 ('appearance', 114),
 ('feeling.', 113),
 ('struggle', 110),
 ('ahh', 109),
 ('struggled', 108),
 ('smelled', 108),
 ('applies', 107),
 ('flawless', 106),
 ('itll', 105),
 ('suggest', 104),
 ('irritating', 104),
 ('dull', 103),
 ('yttp', 102),
 ('better!', 101),
 ('look.', 99),
 ('effects', 98),
 ('applicator', 97),
 ('naturally', 96),
 ('followed', 96),
 ('redness.', 95),
 ('biggest', 93),
 ('stubborn', 93),
 ('sweet', 92),
 ('terrible', 90),
 ('across', 90),
 ('fully', 90),
 ('reallyy', 89),
 ('smoothly', 89),
 ('sunscreens', 89),
 ('affordable', 88),
 ('good!!', 87),
 ('gloss', 87),
 ('addition', 85),
 ('scarring', 84),
 ('setting', 83),
 ('stripping', 83),
 ('disappointed.', 82),
 ('dennis', 82),
 ('weeks,', 80),
 ('too,', 80),
 ('irritation.', 79),
 ('excess', 79),
 ('puffy', 79),
 ('reapply', 77),
 ('better,', 77),
 ('barrier', 77),
 ('amazingg', 76),
 ('initially', 76),
 ('originally', 76),
 ('essential', 76),
 ('smelling', 76),
 ('bottle,', 76),
 ('veryy', 75),
 ('beautifully', 75),
 ('appreciate', 75),
 ('ohh', 74),
 ('tanning', 74),
 ('obsessed.', 73),
 ('pills', 72),
 ('butter', 72),
 ('butt', 71),
 ('currently', 70),
 ('running', 70),
 ('happened', 69),
 ('following', 68),
 ('supple.', 68),
 ('cheeks.', 67),
 ('sunscreen,', 67),
 ('immediately.', 67),
 ('application,', 66),
 ('smoother,', 66),
 ('massage', 66),
 ('impressed.', 66),
 ('allergic', 65),
 ('effective.', 65),
 ('afford', 65),
 ('mirror', 65),
 ('feel.', 65),
 ('difficult', 65),
 ('happy.', 65),
 ('puffiness', 65),
 ('summer.', 64),
 ('worry', 64),
 ('buut', 64),
 ('wanna', 63),
 ('happen', 63),
 ('apply.', 63),
 ('squeeze', 63),
 ('stripped', 63),
 ('accutane', 62),
 ('pill', 61),
 ('difference!', 61),
 ('struggling', 61),
 ('peels', 61),
 ('free.', 61),
 ('stress', 61),
 ('annoying', 61),
 ('bottles', 61),
 ('smoothed', 60),
 ('penny.', 60),
 ('mess', 60),
 ('occasional', 59),
 ('superr', 59),
 ('smoother.', 59),
 ('adds', 59),
 ('pass', 59),
 ('dropper', 58),
 ('cheek', 57),
 ('process', 57),
 ('generally', 57),
 ('deeply', 57),
 ('fall', 57),
 ('yall', 57),
 ('amazingg.', 57),
 ('dollars', 56),
 ('bottom', 56),
 ('issues.', 56),
 ('appear', 56),
 ('obsessed!', 56),
 ('berry', 56),
 ('looved', 55),
 ('employee', 54),
 ('boost', 54),
 ('sk-ii', 53),
 ('beginning', 53),
 ('tanner', 53),
 ('lovee.', 53),
 ('jelly', 53),
 ('glossy', 53),
 ('omgg', 52),
 ('stuff,', 52),
 ('rubbed', 51),
 ('micellar', 51),
 ('gotta', 51),
 ('free,', 51),
 ('cc', 51),
 ('suggested', 50),
 ('verry', 50),
 ('hooked', 50),
 ('sitting', 50),
 ('supergoop', 50),
 ('telling', 50),
 ('agree', 49),
 ('loovvee', 49),
 ('middle', 49),
 ('excellent', 48),
 ('eventually', 48),
 ('smoothing', 48),
 ('current', 48),
 ('yummy', 48),
 ('itt', 47),
 ('loovee', 47),
 ('smaller.', 47),
 ('oiliness', 47),
 ('specifically', 47),
 ('hopefully', 47),
 ('huuge', 47),
 ('ball', 47),
 ('till', 46),
 ('off!', 46),
 ('smell!', 46),
 ('amazingg!', 46),
 ('penny', 45),
 ('whipped', 45),
 ('umm', 45),
 ('strawberry', 45),
 ('biossance', 45),
 ('lovee!', 44),
 ('smaller,', 44),
 ('apply,', 44),
 ('irritated.', 44),
 ('applying.', 44),
 ('mess.', 44),
 ('recommend!!', 44),
 ('dryness.', 44),
 ('needed.', 44),
 ('yellow', 44),
 ('disappoint.', 44),
 ('feet', 43),
 ('reeally', 43),
 ('(especially', 43),
 ('teeth', 43),
 ('school', 42),
 ('odd', 42),
 ('popped', 42),
 ('missing', 41),
 ('good!!!', 41),
 ('smallest', 41),
 ('carry', 41),
 ('impressed!', 41),
 ('buttery', 41),
 ('smooths', 40),
 ('call', 40),
 ('happy!', 40),
 ('hooked.', 40),
 ('scoop', 40),
 ('look,', 40),
 ('vanilla', 40),
 ('suffering', 39),
 ('soothing.', 39),
 ('soothes', 39),
 ('soon.', 39),
 ('effect.', 39),
 ('sunscreen!', 39),
 ('called', 39),
 ('collection', 39),
 ('dramatically', 38),
 ('supple,', 38),
 ('really,', 38),
 ('recommending', 38),
 ('woww', 38),
 ('miss', 38),
 ('sheet', 37),
 ('penny!', 37),
 ('happens', 37),
 ('cheeks,', 37),
 ('planning', 37),
 ('reaally', 37),
 ('possible', 37),
 ('irritates', 37),
 ('dryness,', 36),
 ('suffered', 36),
 ('bottle!', 36),
 ('effective,', 36),
 ('stuff!!', 36),
 ('thankfully', 36),
 ('issue.', 36),
 ('yess', 35),
 ('aand', 35),
 ('messy', 35),
 ('peel.', 35),
 ('practically', 35),
 ('willing', 35),
 ('inflammation', 35),
 ('massaging', 35),
 ('korres', 35),
 ('mm', 34),
 ('goodness', 34),
 ('disappear', 34),
 ('gross.', 34),
 ('deff', 34),
 ('pillow', 34),
 ('soothe', 34),
 ('cells', 34),
 ('looking.', 34),
 ('hmm', 34),
 ('difference,', 33),
 ('little.', 33),
 ('disappeared', 33),
 ('million', 33),
 ('different.', 33),
 ('loove!', 33),
 ('lovvee', 33),
 ('flawless.', 33),
 ('tree', 33),
 ('pull', 33),
 ('correcting', 33),
 ('feeling,', 33),
 ('ooh', 33),
 ('additional', 33),
 ('noo', 32),
 ('thrilled', 32),
 ('applying,', 32),
 ('upper', 32),
 ('bummed', 32),
 ('comment', 32),
 ('sheen', 32),
 ('happier', 32),
 ('steep', 31),
 ('summer,', 31),
 ('brightness', 31),
 ('needing', 31),
 ('combooily', 31),
 ('clogging', 31),
 ('terrible.', 31),
 ('effort', 31),
 ('specially', 31),
 ('lovve', 30),
 ('sorry', 30),
 ('effects.', 30),
 ('fee', 30),
 ('filled', 30),
 ('recommend,', 30),
 ('waayy', 30),
 ('teeny', 30),
 ('collagen', 30),
 ('impossible', 29),
 ('ehh', 29),
 ('loove.', 29),
 ('opportunity', 29),
 ('drastically', 29),
 ('weekly', 29),
 ('looves', 29),
 ('appears', 29),
 ('wonderfully', 29),
 ('gloss.', 29),
 ('disappoint!', 29),
 ('grabbed', 28),
 ('sheer', 28),
 ('thee', 28),
 ('reccomend', 28),
 ('necessarily', 28),
 ('loonngg', 28),
 ('sleep.', 28),
 ('refill', 28),
 ('superfood', 28),
 ('opposed', 28),
 ('bathroom', 28),
 ('hugee', 28),
 ('letting', 28),
 ('freeproduct', 28),
 ('feels.', 27),
 ('irritation,', 27),
 ('commented', 27),
 ('affordable.', 27),
 ('regardless', 27),
 ('andd', 27),
 ('looks.', 27),
 ('applied.', 27),
 ('goop', 27),
 ('pulling', 27),
 ('wallet', 27),
 ('popping', 26),
 ('disappointing', 26),
 ('affect', 26),
 ('disappears', 26),
 ('yess!', 26),
 ('week!', 26),
 ('press', 26),
 ('balls', 26),
 ('missed', 26),
 ('somerville', 26),
 ('ii', 26),
 ('sulwhasoo', 26),
 ('mattifying', 25),
 ('accidentally', 25),
 ('occasionally', 25),
 ('soo,', 25),
 ('lott', 25),
 ('soothed', 25),
 ('dryy', 25),
 ('professional', 25)]

Create correction mapping.

correct_map_1 = {
    r"\bsoo\b": "so",
    r"\bsooo\b": "so",
    r"\bsoooo\b": "so",
    r"\bsooooo\b": "so",

    r"\blovee\b": "love",
    r"\bloove\b": "love",
    r"\bloong\b": "long",
    r"\bwayy\b": "way",
    r"\bgonna\b": "going to",
    r"\bwaay\b": "way",
    r"\bsuuper\b": "super",
    r"\byttp\b": "youth to the people",
    r"\byoull\b": "you will",
    r"\blongg\b": "long",
    r"\bnighttime\b": "night time",    
    r"\bitll\b": "it will",
    r"\breallyy\b": "really", 
    r"\bwanna\b": "want to", 
    r"\bohh\b": " ", 
    r"\bbuut\b": "but", 
    r"\byall\b": "you all",
    r"\bloovvee\b": "love",    
    r"\blovee\.": "love",      
    r"\bamazingg\b": "amazing",       
    r"\bwoww\b": "wow",     
    r"\bsuperr\b": "super",
    r"\baand\b": "and",
    r"\bamazingg\.": "amazing", 
    r"\bnoo\b": "no", 
    r"\bgotta\b": "got to", 
    r"\bdeff\b": "definitely", 
    r"\bveryy\b": "very", 
    r"\bandd\b": "and",
    r"\bomgg\b": "oh my god",
    r"\bloove\.": "love",
    r"\bloovee\b": "love",
    r"\blooved\b": "loved",
    r"\byess\b": "yes",
    r"\byall,": "you all",
    r"\bhuuge\b": "huge",
    r"\byess!": "yes",
    r"\blovvee\b": "love",
    r"\bwaayy\b": "way",
    r"\bloove!\b": "love",
    r"\bbutt\b": "but",
    r"\bloonngg\b": "long",
    r"\bitt\.": "it",
    r"\bforeverr\b": "forever",
    r"\bthicc\b": "thick",
    r"\bverry\b": "very",
    r"\bitt\b": "it",
    r"\blonng\b": "long",
    r"\beverr\b": "ever",
    r"\booh\b": " ", 
    r"\blovedd\b": "loved",
    r"\blott\b": "lot",
    r"\bsoo,": "so",
    r"\bmuchh\b": "much",
    r"\byess!!!": "yes",
    r"\bwoow\b": " ",
    r"\blooves\b": "loves",
    r"\byall\.": "you all",
    r"\bitt!": "it",
    r"\breaally\b": "really",
    r"\bamaazing\.": "amazing", 
    r"\bverryy\b": "very",
    r"\bhella\b": "very",
    r"\bamazingg!": "amazing", 
    r"\bloovvee\b": "love",

    # later added
    r"\bahh\b": " ", 
    r"\bumm\b": " ",     
    r"\blovee!": "love",
    r"\byess\b": "yes",
    r"\bhmm\b": " ", 
    r"\bcombooily\b": "combo oily", 
    r"\blovve\b": "love",   
    r"\behh\b": " ",  
    r"\bloonngg\b": "long",    
    r"\bhugee\b": "huge",  
    r"\bfreeproduct\b": "free product",    
    r"\bandd\b": "and",
    r"\byess!": "yes",
    r"\blott\b": "lot",
}

review_df_all_deduplicate['review_text_cleaned'] = (
    review_df_all_deduplicate['review_text_cleaned']
    .replace(correct_map_1,regex=True)
)

4.8 Punctuation

Remove punctuations.

pattern = "[" + re.escape(string.punctuation) + "]" 
review_df_all_deduplicate["review_text_cleaned"] = (
    review_df_all_deduplicate["review_text_cleaned"].str.replace(pattern, " ", regex=True)
)

punc_lis_remain = []

for text in review_df_all_deduplicate['review_text_cleaned']:
    punc_lis_remain.extend(re.findall(r"[^\w\s]", text))

punc_lis_remain = sorted(set(punc_lis_remain), key=lambda ch: ord(ch)) 
print('Characters (only emoji) left in the reviews: ', '\n', punc_lis_remain)

Characters (only emoji) left in the reviews:  
 ['☀', '☁', '★', '☆', '☝', '☬', '☹', '☺', '♀', '♂', '♡', '♥', '♻', '♾', '⚗', '⚜', '⚠', '⚡', '⚪', '⛅', '⛑', '✅', '✈', '✋', '✌', '✓', '✔', '✖', '✨', '❄', '❌', '❕', '❗', '❣', '❤', '➕', '➖', '➡', '➾', '➿', '⬇', '⭐', '�', '🅶', '🅷', '🅼', '🅾', '🆂', '🆃', '🆄', '🆈', '🆓', '🆙', '🌊', '🌙', '🌞', '🌟', '🌫', '🌱', '🌸', '🌹', '🌺', '🌿', '🍂', '🍃', '🍅', '🍉', '🍊', '🍋', '🍌', '🍒', '🍓', '🍞', '🍬', '🍯', '🎀', '🎄', '🎉', '🏆', '🏻', '🏼', '🏽', '🏾', '🐐', '🐘', '🐝', '🐣', '👀', '👁', '👄', '👋', '👌', '👍', '👎', '👏', '👣', '👩', '👶', '💀', '💄', '💅', '💋', '💓', '💔', '💕', '💖', '💗', '💘', '💙', '💚', '💛', '💜', '💞', '💡', '💥', '💦', '💨', '💫', '💯', '📌', '📍', '📦', '🔥', '🔵', '🖤', '😀', '😁', '😂', '😃', '😄', '😅', '😆', '😉', '😊', '😌', '😍', '😏', '😒', '😓', '😔', '😕', '😖', '😘', '😚', '😝', '😞', '😠', '😢', '😩', '😪', '😫', '😬', '😭', '😮', '😲', '😳', '😶', '😻', '🙁', '🙂', '🙃', '🙄', '🙌', '🙏', '🛍', '🛑', '🛒', '🤌', '🤍', '🤎', '🤐', '🤓', '🤔', '🤗', '🤞', '🤡', '🤢', '🤣', '🤤', '🤦', '🤨', '🤩', '🤪', '🤭', '🤮', '🤯', '🤷', '🥑', '🥰', '🥲', '🥳', '🥴', '🥶', '🥹', '🥺', '🦝', '🦲', '🧖', '🧡', '🧴', '🧼', '🧿', '🩸', '🪄', '🫠', '🫣', '🫰', '🫶']

4.9 Random Words

Finally, notice the existence of random and meaningless words that have extremely long length, for example:

review_df_all_deduplicate.loc[80989, 'review_text_cleaned']

'yea i agree yes same mhm jtssngdbkgiiycycitctiheckljjioo'

We applied a heuristic threshold of 23 characters, any individual token exceeding this limit was identified as an anomaly and removed.

def detect_gibberish(w, min_len = 30):    
    if len(w) > min_len:
        return True

    return False

def gibberish_text(text):
    words = text.split()
    return any(detect_gibberish(w) for w in words)


mask = review_df_all_deduplicate['review_text_cleaned'].apply(gibberish_text)

df_process = review_df_all_deduplicate.loc[mask, 'review_text_cleaned']
df_process = pd.DataFrame(df_process)

def drop_long_word(text, threshold = 23):
    emoji = {"✅", "❌"} 
    tokens = text.split()
    kept = []

    for tok in tokens:
        if len(tok) >= threshold:
            if any(e in tok for e in emoji):
                kept.append(tok)

        else:
            kept.append(tok)

    return " ".join(kept)


df_process['review_text_cleaned'] = df_process['review_text_cleaned'].apply(drop_long_word)
review_df_all_deduplicate.loc[df_process.index, 'review_text_cleaned'] = df_process['review_text_cleaned']

This concludes our text cleaning pipeline. While achieving a perfectly clean corpus is impractical given the large volume of unstructured text, this process has substantially reduced noise and standardized the data.

review_df_all_deduplicate['review_text_cleaned'] = review_df_all_deduplicate['review_text_cleaned'].str.replace(r"\s+", " ", regex=True)
review_df_all_deduplicate['review_text_cleaned'] = review_df_all_deduplicate['review_text_cleaned'].str.strip()

# some text entries are empty after cleaning, we remove them
idx_length_0 = review_df_all_deduplicate[review_df_all_deduplicate['review_text_cleaned'].map(len) == 0].index

mask = np.isin(review_df_all_deduplicate.index, idx_length_0)
review_df_all_deduplicate.loc[review_df_all_deduplicate.index.values[~ mask]].to_csv('dataset/review_df_all_deduplicate_processed.csv', index=True)

5. Sentence Embedding

We use Sentence-BERT to generate sentence embeddings. Unlike traditional keyword-based methods (e.g. TF-IDF) or static embeddings (e.g. Fasttext), SBERT captures deep contextual meaning. For instance, it understands that “not good” is the opposite of “good,” whereas keyword-based models might treat them similarly. Therefore it is the most suitable model to precisely convert sentences into vectors given our current computational resources and data.

The implementation details are omitted for brevity. Since this model is a transformer-based model that benefits significantly from GPU acceleration, the full execution script is hosted on Google Colab and can be found here.

We directly import the transformed vectors.

df_sbert = pd.read_csv("dataset/sbert_embeddings_2.csv", index_col=0)

6. Clustering

Given our dataset size of over 800,000 reviews, performing t-SNE to firstly visualize clustering structure of the dataset is computationally prohibitive, we use UMAP instead, since:

1. UMAP is much faster especially on large datasets.

2. UMAP better maintains the relationships between clusters, not just points within clusters, offering a more accurate representation of the high-dimensional data’s overall shape.

6.1 UMAP

We first draw a set of UMAP plots with different combinations of n_neighbors and min_dist to avoid instability caused by parameters. Candidate values are n_neighbors = 15, 30, 45, 60, and min_dist = 0.1.

Notice here we choose the init = ‘pca’ and metric = ‘cosine’.

n_neighbors = 15, min_dist = 0.1.

%time
umap = UMAP(n_neighbors=15, min_dist = 0.1, n_components=2, init='pca', metric='cosine', n_jobs=-1, verbose=True)
mat_umap = umap.fit_transform(df_sbert)

CPU times: total: 0 ns
Wall time: 0 ns
UMAP(angular_rp_forest=True, init='pca', metric='cosine', verbose=True)
Fri Jan  9 22:31:20 2026 Construct fuzzy simplicial set
Fri Jan  9 22:31:21 2026 Finding Nearest Neighbors
Fri Jan  9 22:31:21 2026 Building RP forest with 52 trees
Fri Jan  9 22:32:17 2026 NN descent for 20 iterations
     1  /  20
     2  /  20
     3  /  20
     4  /  20
     5  /  20
     6  /  20
     7  /  20
    Stopping threshold met -- exiting after 7 iterations
Fri Jan  9 22:33:22 2026 Finished Nearest Neighbor Search
Fri Jan  9 22:33:32 2026 Construct embedding

    completed  0  /  200 epochs
    completed  20  /  200 epochs
    completed  40  /  200 epochs
    completed  60  /  200 epochs
    completed  80  /  200 epochs
    completed  100  /  200 epochs
    completed  120  /  200 epochs
    completed  140  /  200 epochs
    completed  160  /  200 epochs
    completed  180  /  200 epochs
Fri Jan  9 22:35:37 2026 Finished embedding

n_neighbors = 30, min_dist = 0.1.

%time
umap = UMAP(n_neighbors=30, min_dist = 0.1, n_components=2, init='pca', metric='cosine', n_jobs=-1, verbose=True)
mat_umap_1 = umap.fit_transform(df_sbert)

CPU times: total: 0 ns
Wall time: 0 ns
UMAP(angular_rp_forest=True, init='pca', metric='cosine', n_neighbors=30, verbose=True)
Fri Jan  9 22:36:35 2026 Construct fuzzy simplicial set
Fri Jan  9 22:36:37 2026 Finding Nearest Neighbors
Fri Jan  9 22:36:37 2026 Building RP forest with 52 trees
Fri Jan  9 22:37:27 2026 NN descent for 20 iterations
     1  /  20
     2  /  20
     3  /  20
     4  /  20
     5  /  20
     6  /  20
    Stopping threshold met -- exiting after 6 iterations
Fri Jan  9 22:40:12 2026 Finished Nearest Neighbor Search
Fri Jan  9 22:40:22 2026 Construct embedding

    completed  0  /  200 epochs
    completed  20  /  200 epochs
    completed  40  /  200 epochs
    completed  60  /  200 epochs
    completed  80  /  200 epochs
    completed  100  /  200 epochs
    completed  120  /  200 epochs
    completed  140  /  200 epochs
    completed  160  /  200 epochs
    completed  180  /  200 epochs
Fri Jan  9 22:43:32 2026 Finished embedding

n_neighbors = 45, min_dist = 0.1.

%time
umap = UMAP(n_neighbors=45, min_dist = 0.1, n_components=2, init='pca', metric='cosine', n_jobs=-1, verbose=True)
mat_umap_2 = umap.fit_transform(df_sbert)

CPU times: total: 0 ns
Wall time: 0 ns
UMAP(angular_rp_forest=True, init='pca', metric='cosine', n_neighbors=45, verbose=True)
Tue Jan 13 00:27:07 2026 Construct fuzzy simplicial set
Tue Jan 13 00:27:08 2026 Finding Nearest Neighbors
Tue Jan 13 00:27:09 2026 Building RP forest with 52 trees
Tue Jan 13 00:28:15 2026 NN descent for 20 iterations
     1  /  20
     2  /  20
     3  /  20
     4  /  20
     5  /  20
    Stopping threshold met -- exiting after 5 iterations
Tue Jan 13 00:30:55 2026 Finished Nearest Neighbor Search
Tue Jan 13 00:31:10 2026 Construct embedding

    completed  0  /  200 epochs
    completed  20  /  200 epochs
    completed  40  /  200 epochs
    completed  60  /  200 epochs
    completed  80  /  200 epochs
    completed  100  /  200 epochs
    completed  120  /  200 epochs
    completed  140  /  200 epochs
    completed  160  /  200 epochs
    completed  180  /  200 epochs
Tue Jan 13 00:33:56 2026 Finished embedding

n_neighbors = 60, min_dist = 0.1.

%time
umap = UMAP(n_neighbors=60, min_dist = 0.1, n_components=2, init='pca', metric='cosine', n_jobs=-1, verbose=True)
mat_umap_3 = umap.fit_transform(df_sbert)

CPU times: total: 0 ns
Wall time: 0 ns
UMAP(angular_rp_forest=True, init='pca', metric='cosine', n_neighbors=60, verbose=True)
Fri Jan  9 22:53:26 2026 Construct fuzzy simplicial set
Fri Jan  9 22:53:27 2026 Finding Nearest Neighbors
Fri Jan  9 22:53:27 2026 Building RP forest with 52 trees
Fri Jan  9 22:54:10 2026 NN descent for 20 iterations
     1  /  20
     2  /  20
     3  /  20
     4  /  20
     5  /  20
    Stopping threshold met -- exiting after 5 iterations
Fri Jan  9 23:00:28 2026 Finished Nearest Neighbor Search
Fri Jan  9 23:00:46 2026 Construct embedding

    completed  0  /  200 epochs
    completed  20  /  200 epochs
    completed  40  /  200 epochs
    completed  60  /  200 epochs
    completed  80  /  200 epochs
    completed  100  /  200 epochs
    completed  120  /  200 epochs
    completed  140  /  200 epochs
    completed  160  /  200 epochs
    completed  180  /  200 epochs
Fri Jan  9 23:05:10 2026 Finished embedding

fig, axes = plt.subplots(2,2, figsize = (16, 12))

df_umap = pd.DataFrame(mat_umap)
sns.scatterplot(x=df_umap.iloc[:,0], y=df_umap.iloc[:,1], alpha = 0.2, ax = axes[0,0])
axes[0,0].set_xlabel("x_projected")
axes[0,0].set_ylabel("y_projected")
axes[0,0].set_title('n_neighbors: 15, min_dist = 0.1')


df_umap = pd.DataFrame(mat_umap_1)
sns.scatterplot(x=df_umap.iloc[:,0], y=df_umap.iloc[:,1], alpha = 0.2, ax = axes[0,1])
axes[0,1].set_xlabel("x_projected")
axes[0,1].set_ylabel("y_projected")
axes[0,1].set_title('n_neighbors: 30, min_dist = 0.1')

df_umap = pd.DataFrame(mat_umap_2)
sns.scatterplot(x=df_umap.iloc[:,0], y=df_umap.iloc[:,1], alpha = 0.2, ax = axes[1,0])
axes[1,0].set_xlabel("x_projected")
axes[1,0].set_ylabel("y_projected")
axes[1,0].set_title('n_neighbors: 45, min_dist = 0.1')

df_umap = pd.DataFrame(mat_umap_3)
sns.scatterplot(x=df_umap.iloc[:,0], y=df_umap.iloc[:,1], alpha = 0.2, ax = axes[1,1])
axes[1,1].set_xlabel("x_projected")
axes[1,1].set_ylabel("y_projected")
axes[1,1].set_title('n_neighbors: 60, min_dist = 0.1')

fig.suptitle('UMAP', fontsize = 24)

plt.tight_layout()
plt.show()

Select n_neighbors = 45, min_dist = 0.1 as our representative umap plot.

fig, ax = plt.subplots()
df_umap = pd.DataFrame(mat_umap_2)
sns.scatterplot(x=df_umap.iloc[:,0], y=df_umap.iloc[:,1], alpha = 0.2)
ax.set_xlabel("x_projected")
ax.set_ylabel("y_projected")
ax.set_title('n_neighbors: 45, min_dist = 0.1')
plt.show()

We see a main larger cluster surronded by 3 dense smaller clusters, they are not well-separated. Two of the smaller clusters are connected within bridge points, and additionally considering the shape of the clusters, they seem not convex, algorithms like k-means and hierarchical clustering are not suitable to our case.

6.1.1 Reduced Dimension

Before implementing the clustering algorithm, dimensionality reduction step is needed, as the sbert embedded vectors in 384-dimension are suffering ‘curse of dimensionality’, the distance between any 2 points converges and becomes nearly equidistant in the high dimensional space. By reducing dimensions to a latent space, UMAP concentrates the variance and re-establishes a meaningful distance metric, making it easier for algorithms to identify distinct boundaries.

We use the trustworthiness to evaluate optimal reducing dimension. The parameter n_neighbors needs to be tuned for the calculation of cluster trustworthiness, four values (15, 30, 50, 100) of n_neighbors will be tested, a small value of n_neighbors tends to preserve local clustering structure, while a larger value tends to preserve the global clustering structure, avoiding overffing.

df_sbert['rating'] = review_df_all_deduplicate.loc[df_sbert.index, 'rating']

sample_ratio = 0.023

_, df_sbert_sample = train_test_split(df_sbert, test_size=sample_ratio, stratify=df_sbert['rating'], 
                                random_state=42)

def umap_dimension_select(n_neighbors, dim_start=20, dim_end=51, df = df_sbert_sample.iloc[:, : -1]):
    """
    pandas df
    """

    dim_num_list = range(dim_start, dim_end)
    trustworthiness_list=[]

    for k in dim_num_list:
        
        mat = df.to_numpy()
        
        umap = UMAP(n_neighbors=45, min_dist = 0.1, n_components=k, init='pca', metric='cosine', n_jobs=-1, verbose=False)
        mat_reduced_sample = umap.fit_transform(mat)
        
        trust_score = trustworthiness(mat, mat_reduced_sample, n_neighbors=n_neighbors)
        
        trustworthiness_list.append(trust_score)

    return trustworthiness_list

fig, axes = plt.subplots(2,2,figsize = (12, 8))
neighbor_lis = [15, 30, 50, 100]

dim_num_list = range(20, 51)

for neighbor_num in neighbor_lis:
     
     if neighbor_num == 15:
        
        trustworthiness_list = umap_dimension_select(n_neighbors=neighbor_num)
        
        print(f"neighbor number {neighbor_num} completed.")
        
        axes[0,0].plot(dim_num_list, trustworthiness_list)
        axes[0,0].set_xlabel('UMAP dimension')
        axes[0,0].set_xticks(list(dim_num_list))
        axes[0,0].tick_params(axis='x', rotation=45, labelsize=10)
        axes[0,0].set_ylabel('Trustworthiness score')
        axes[0,0].set_title(f'Trustworthiness Score Plot (neighbor_num: {neighbor_num})')
        axes[0,0].grid(True, alpha=0.3)

     elif neighbor_num == 30:
        
        trustworthiness_list = umap_dimension_select(n_neighbors=neighbor_num)
        
        print(f"neighbor number {neighbor_num} completed.")
        
        axes[0,1].plot(dim_num_list, trustworthiness_list)
        axes[0,1].set_xlabel('UMAP dimension')
        axes[0,1].set_xticks(list(dim_num_list))
        axes[0,1].tick_params(axis='x', rotation=45, labelsize=10)
        axes[0,1].set_ylabel('Trustworthiness score')
        axes[0,1].set_title(f'Trustworthiness Score Plot (neighbor_num: {neighbor_num})')
        axes[0,1].grid(True, alpha=0.3)

     elif neighbor_num == 50:
        
        trustworthiness_list = umap_dimension_select(n_neighbors=neighbor_num)
        
        print(f"neighbor number {neighbor_num} completed.")
        
        axes[1,0].plot(dim_num_list, trustworthiness_list)
        axes[1,0].set_xlabel('UMAP dimension')
        axes[1,0].set_xticks(list(dim_num_list))
        axes[1,0].tick_params(axis='x', rotation=45, labelsize=10)
        axes[1,0].set_ylabel('Trustworthiness score')
        axes[1,0].set_title(f'Trustworthiness Score Plot (neighbor_num: {neighbor_num})')
        axes[1,0].grid(True, alpha=0.3)

     else:
          
        trustworthiness_list = umap_dimension_select(n_neighbors=neighbor_num)
        
        print(f"neighbor number {neighbor_num} completed.")
        
        axes[1,1].plot(dim_num_list, trustworthiness_list)
        axes[1,1].set_xlabel('UMAP dimension')
        axes[1,1].set_xticks(list(dim_num_list))
        axes[1,1].tick_params(axis='x', rotation=45, labelsize=10)
        axes[1,1].set_ylabel('Trustworthiness score')
        axes[1,1].set_title(f'Trustworthiness Score Plot (neighbor_num: {neighbor_num})')
        axes[1,1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

neighbor number 15 completed.
neighbor number 30 completed.
neighbor number 50 completed.
neighbor number 100 completed.

The scale of the score is significantly lower than the other 3 n_neighbors when its value is 100, it included too many unnecessary neighbors that resulted in a high penalty. And in general, the other 3 plots agree reducing the original dimension to 50 is a solid choice, even though it did not score the best when n_neighbors are 15 and 30.

umap = UMAP(n_neighbors=45, min_dist = 0.1, n_components=50, init='pca', metric='cosine', n_jobs=-1, verbose=True)
mat_umap_cluster = umap.fit_transform(df_sbert)

UMAP(angular_rp_forest=True, init='pca', metric='cosine', n_components=50, n_neighbors=45, verbose=True)
Tue Jan 13 00:34:17 2026 Construct fuzzy simplicial set
Tue Jan 13 00:34:18 2026 Finding Nearest Neighbors
Tue Jan 13 00:34:19 2026 Building RP forest with 52 trees
Tue Jan 13 00:35:01 2026 NN descent for 20 iterations
     1  /  20
     2  /  20
     3  /  20
     4  /  20
     5  /  20
    Stopping threshold met -- exiting after 5 iterations
Tue Jan 13 00:37:46 2026 Finished Nearest Neighbor Search
Tue Jan 13 00:38:02 2026 Construct embedding

    completed  0  /  200 epochs
    completed  20  /  200 epochs
    completed  40  /  200 epochs
    completed  60  /  200 epochs
    completed  80  /  200 epochs
    completed  100  /  200 epochs
    completed  120  /  200 epochs
    completed  140  /  200 epochs
    completed  160  /  200 epochs
    completed  180  /  200 epochs
Tue Jan 13 00:44:16 2026 Finished embedding

umap = UMAP(n_neighbors=45, min_dist = 0.1, n_components=50, init='pca', metric='cosine', n_jobs=-1, verbose=True)
mat_umap_cluster_sample = umap.fit_transform(df_sbert_sample.iloc[:, : -1])

UMAP(angular_rp_forest=True, init='pca', metric='cosine', n_components=50, n_neighbors=45, verbose=True)
Tue Jan 13 00:44:24 2026 Construct fuzzy simplicial set
Tue Jan 13 00:44:24 2026 Finding Nearest Neighbors
Tue Jan 13 00:44:24 2026 Building RP forest with 12 trees
Tue Jan 13 00:44:24 2026 NN descent for 14 iterations
     1  /  14
     2  /  14
     3  /  14
     4  /  14
    Stopping threshold met -- exiting after 4 iterations
Tue Jan 13 00:44:27 2026 Finished Nearest Neighbor Search
Tue Jan 13 00:44:27 2026 Construct embedding

    completed  0  /  200 epochs
    completed  20  /  200 epochs
    completed  40  /  200 epochs
    completed  60  /  200 epochs
    completed  80  /  200 epochs
    completed  100  /  200 epochs
    completed  120  /  200 epochs
    completed  140  /  200 epochs
    completed  160  /  200 epochs
    completed  180  /  200 epochs
Tue Jan 13 00:44:35 2026 Finished embedding

6.2 HDBSCAN

def hdbscan_test(min_cluster_size, min_samples = 10, df = mat_umap_cluster_sample):
    hdb = hdbscan.HDBSCAN(min_cluster_size= min_cluster_size, min_samples=min_samples)
    hdb.fit(df)

    print('min_cluster_size:', min_cluster_size, ',', 'min_samples:', min_samples)
    print('number of cluster: ', len(np.unique(hdb.labels_)))
    print('noise proportion: ', len(hdb.labels_[hdb.labels_ == -1]) / len(hdb.labels_))  

hdbscan_test(min_cluster_size=50)
print(' ')
hdbscan_test(min_cluster_size=80)
print(' ')
hdbscan_test(min_cluster_size=300)

min_cluster_size: 50 , min_samples: 10
number of cluster:  43
noise proportion:  0.4299823978095052
 
min_cluster_size: 80 , min_samples: 10
number of cluster:  33
noise proportion:  0.3624095442988461
 
min_cluster_size: 300 , min_samples: 10
number of cluster:  5
noise proportion:  0.007138666145120281

min_cluster_size_range = np.arange(50, 301, 1) 
min_samples = 10

silhouette_scores = []
calinski_scores = []
davies_scores = []
    
num_clusters_list = []
noise_proportion_list = []

for min_cluster in min_cluster_size_range:


    hdb = hdbscan.HDBSCAN(min_cluster_size= min_cluster, min_samples=min_samples)
    hdb.fit(mat_umap_cluster_sample)

    pred_label = hdb.labels_

    n_clusters = len(set(pred_label)) - (1 if -1 in pred_label else 0)
    n_noise = list(pred_label).count(-1)

    num_clusters_list.append(n_clusters)
    noise_proportion_list.append(n_noise / len(mat_umap_cluster_sample))

    if n_clusters > 1 and n_noise < len(mat_umap_cluster_sample):
        mask = pred_label != -1
        silhou_score = silhouette_score(mat_umap_cluster_sample[mask], pred_label[mask])
        silhouette_scores.append(silhou_score)

        calinski_score = calinski_harabasz_score(mat_umap_cluster_sample[mask], pred_label[mask])
        calinski_scores.append(calinski_score)

        davies_score = davies_bouldin_score(mat_umap_cluster_sample[mask], pred_label[mask])
        davies_scores.append(davies_score)

    else:
        silhouette_scores.append(0)
        calinski_scores.append(0)
        davies_scores.append(0)
    
    print(f"min_cluster: {min_cluster} completed.")

min_cluster: 50 completed.
min_cluster: 51 completed.
min_cluster: 52 completed.
min_cluster: 53 completed.
min_cluster: 54 completed.
min_cluster: 55 completed.
min_cluster: 56 completed.
min_cluster: 57 completed.
min_cluster: 58 completed.
min_cluster: 59 completed.
min_cluster: 60 completed.
min_cluster: 61 completed.
min_cluster: 62 completed.
min_cluster: 63 completed.
min_cluster: 64 completed.
min_cluster: 65 completed.
min_cluster: 66 completed.
min_cluster: 67 completed.
min_cluster: 68 completed.
min_cluster: 69 completed.
min_cluster: 70 completed.
min_cluster: 71 completed.
min_cluster: 72 completed.
min_cluster: 73 completed.
min_cluster: 74 completed.
min_cluster: 75 completed.
min_cluster: 76 completed.
min_cluster: 77 completed.
min_cluster: 78 completed.
min_cluster: 79 completed.
min_cluster: 80 completed.
min_cluster: 81 completed.
min_cluster: 82 completed.
min_cluster: 83 completed.
min_cluster: 84 completed.
min_cluster: 85 completed.
min_cluster: 86 completed.
min_cluster: 87 completed.
min_cluster: 88 completed.
min_cluster: 89 completed.
min_cluster: 90 completed.
min_cluster: 91 completed.
min_cluster: 92 completed.
min_cluster: 93 completed.
min_cluster: 94 completed.
min_cluster: 95 completed.
min_cluster: 96 completed.
min_cluster: 97 completed.
min_cluster: 98 completed.
min_cluster: 99 completed.
min_cluster: 100 completed.
min_cluster: 101 completed.
min_cluster: 102 completed.
min_cluster: 103 completed.
min_cluster: 104 completed.
min_cluster: 105 completed.
min_cluster: 106 completed.
min_cluster: 107 completed.
min_cluster: 108 completed.
min_cluster: 109 completed.
min_cluster: 110 completed.
min_cluster: 111 completed.
min_cluster: 112 completed.
min_cluster: 113 completed.
min_cluster: 114 completed.
min_cluster: 115 completed.
min_cluster: 116 completed.
min_cluster: 117 completed.
min_cluster: 118 completed.
min_cluster: 119 completed.
min_cluster: 120 completed.
min_cluster: 121 completed.
min_cluster: 122 completed.
min_cluster: 123 completed.
min_cluster: 124 completed.
min_cluster: 125 completed.
min_cluster: 126 completed.
min_cluster: 127 completed.
min_cluster: 128 completed.
min_cluster: 129 completed.
min_cluster: 130 completed.
min_cluster: 131 completed.
min_cluster: 132 completed.
min_cluster: 133 completed.
min_cluster: 134 completed.
min_cluster: 135 completed.
min_cluster: 136 completed.
min_cluster: 137 completed.
min_cluster: 138 completed.
min_cluster: 139 completed.
min_cluster: 140 completed.
min_cluster: 141 completed.
min_cluster: 142 completed.
min_cluster: 143 completed.
min_cluster: 144 completed.
min_cluster: 145 completed.
min_cluster: 146 completed.
min_cluster: 147 completed.
min_cluster: 148 completed.
min_cluster: 149 completed.
min_cluster: 150 completed.
min_cluster: 151 completed.
min_cluster: 152 completed.
min_cluster: 153 completed.
min_cluster: 154 completed.
min_cluster: 155 completed.
min_cluster: 156 completed.
min_cluster: 157 completed.
min_cluster: 158 completed.
min_cluster: 159 completed.
min_cluster: 160 completed.
min_cluster: 161 completed.
min_cluster: 162 completed.
min_cluster: 163 completed.
min_cluster: 164 completed.
min_cluster: 165 completed.
min_cluster: 166 completed.
min_cluster: 167 completed.
min_cluster: 168 completed.
min_cluster: 169 completed.
min_cluster: 170 completed.
min_cluster: 171 completed.
min_cluster: 172 completed.
min_cluster: 173 completed.
min_cluster: 174 completed.
min_cluster: 175 completed.
min_cluster: 176 completed.
min_cluster: 177 completed.
min_cluster: 178 completed.
min_cluster: 179 completed.
min_cluster: 180 completed.
min_cluster: 181 completed.
min_cluster: 182 completed.
min_cluster: 183 completed.
min_cluster: 184 completed.
min_cluster: 185 completed.
min_cluster: 186 completed.
min_cluster: 187 completed.
min_cluster: 188 completed.
min_cluster: 189 completed.
min_cluster: 190 completed.
min_cluster: 191 completed.
min_cluster: 192 completed.
min_cluster: 193 completed.
min_cluster: 194 completed.
min_cluster: 195 completed.
min_cluster: 196 completed.
min_cluster: 197 completed.
min_cluster: 198 completed.
min_cluster: 199 completed.
min_cluster: 200 completed.
min_cluster: 201 completed.
min_cluster: 202 completed.
min_cluster: 203 completed.
min_cluster: 204 completed.
min_cluster: 205 completed.
min_cluster: 206 completed.
min_cluster: 207 completed.
min_cluster: 208 completed.
min_cluster: 209 completed.
min_cluster: 210 completed.
min_cluster: 211 completed.
min_cluster: 212 completed.
min_cluster: 213 completed.
min_cluster: 214 completed.
min_cluster: 215 completed.
min_cluster: 216 completed.
min_cluster: 217 completed.
min_cluster: 218 completed.
min_cluster: 219 completed.
min_cluster: 220 completed.
min_cluster: 221 completed.
min_cluster: 222 completed.
min_cluster: 223 completed.
min_cluster: 224 completed.
min_cluster: 225 completed.
min_cluster: 226 completed.
min_cluster: 227 completed.
min_cluster: 228 completed.
min_cluster: 229 completed.
min_cluster: 230 completed.
min_cluster: 231 completed.
min_cluster: 232 completed.
min_cluster: 233 completed.
min_cluster: 234 completed.
min_cluster: 235 completed.
min_cluster: 236 completed.
min_cluster: 237 completed.
min_cluster: 238 completed.
min_cluster: 239 completed.
min_cluster: 240 completed.
min_cluster: 241 completed.
min_cluster: 242 completed.
min_cluster: 243 completed.
min_cluster: 244 completed.
min_cluster: 245 completed.
min_cluster: 246 completed.
min_cluster: 247 completed.
min_cluster: 248 completed.
min_cluster: 249 completed.
min_cluster: 250 completed.
min_cluster: 251 completed.
min_cluster: 252 completed.
min_cluster: 253 completed.
min_cluster: 254 completed.
min_cluster: 255 completed.
min_cluster: 256 completed.
min_cluster: 257 completed.
min_cluster: 258 completed.
min_cluster: 259 completed.
min_cluster: 260 completed.
min_cluster: 261 completed.
min_cluster: 262 completed.
min_cluster: 263 completed.
min_cluster: 264 completed.
min_cluster: 265 completed.
min_cluster: 266 completed.
min_cluster: 267 completed.
min_cluster: 268 completed.
min_cluster: 269 completed.
min_cluster: 270 completed.
min_cluster: 271 completed.
min_cluster: 272 completed.
min_cluster: 273 completed.
min_cluster: 274 completed.
min_cluster: 275 completed.
min_cluster: 276 completed.
min_cluster: 277 completed.
min_cluster: 278 completed.
min_cluster: 279 completed.
min_cluster: 280 completed.
min_cluster: 281 completed.
min_cluster: 282 completed.
min_cluster: 283 completed.
min_cluster: 284 completed.
min_cluster: 285 completed.
min_cluster: 286 completed.
min_cluster: 287 completed.
min_cluster: 288 completed.
min_cluster: 289 completed.
min_cluster: 290 completed.
min_cluster: 291 completed.
min_cluster: 292 completed.
min_cluster: 293 completed.
min_cluster: 294 completed.
min_cluster: 295 completed.
min_cluster: 296 completed.
min_cluster: 297 completed.
min_cluster: 298 completed.
min_cluster: 299 completed.

fig, axes = plt.subplots(2, 3, figsize = (25, 10))

axes[0,0].plot(min_cluster_size_range,silhouette_scores)
axes[0,0].set_xlabel('min_cluster size')
axes[0,0].set_ylabel('silhouette score')
axes[0,0].set_title('silhouette score plot')
axes[0,0].set_xticks(list(range(50, 300, 5)))
axes[0,0].grid(True, alpha = 0.1)
axes[0,0].tick_params(axis='x', rotation = 45, labelsize = 8)
axes[0,0].axvline(x = 194, color = 'red', linestyle = '--')

axes[0,1].plot(min_cluster_size_range,calinski_scores)
axes[0,1].set_xlabel('min_cluster size')
axes[0,1].set_ylabel('calinski score')
axes[0,1].set_title('calinski score plot')
axes[0,1].set_xticks(list(range(50, 300, 5)))
axes[0,1].grid(True, alpha = 0.1)
axes[0,1].tick_params(axis='x', rotation = 45, labelsize = 8)

axes[0,2].plot(min_cluster_size_range,davies_scores)
axes[0,2].set_xlabel('min_cluster size')
axes[0,2].set_ylabel('davies score')
axes[0,2].set_title('davies score plot')
axes[0,2].set_xticks(list(range(50, 300, 5)))
axes[0,2].grid(True, alpha = 0.1)
axes[0,2].tick_params(axis='x', rotation = 45, labelsize = 8)

axes[1,0].plot(min_cluster_size_range,num_clusters_list)
axes[1,0].set_xlabel('min_cluster size')
axes[1,0].set_ylabel('cluster number')
axes[1,0].set_title('cluster number plot')
axes[1,0].set_xticks(list(range(50, 300, 5)))
axes[1,0].grid(True, alpha = 0.1)
axes[1,0].tick_params(axis='x', rotation = 45, labelsize = 8)

axes[1,1].plot(min_cluster_size_range, noise_proportion_list)
axes[1,1].set_xlabel('min_cluster size')
axes[1,1].set_ylabel('noise proportion')
axes[1,1].set_title('noise plot')
axes[1,1].set_xticks(list(range(50, 300, 5)))
axes[1,1].grid(True, alpha = 0.1)
axes[1,1].tick_params(axis='x', rotation = 45, labelsize = 8)

axes[1,2].axis('off')
plt.tight_layout()
plt.show()

Implement on the full dataset.

hdb = hdbscan.HDBSCAN(min_cluster_size= 194, min_samples= 10)
hdb.fit(mat_umap_cluster)

6.1 PCA + K-Means

6.1.1 PCA

mat = df_sbert.to_numpy()

mat_normalized = normalize(mat, norm='l2')
pca = PCA(n_components=200, random_state= 42).fit(mat_normalized)
cumvar = np.cumsum(pca.explained_variance_ratio_)

plt.figure(figsize=(6,4))
plt.plot(range(1, len(cumvar)+1), cumvar, marker="o")
plt.axhline(0.9, color="r", ls="--", lw=1) 
plt.xlabel("Number of components")
plt.ylabel("Cumulative explained variance")
plt.ylim(0, 1.01)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print('principal component number containing 90% variance: ',np.argmin(cumvar <= 0.9) + 1)

principal component number containing 90% variance:  152

pca = PCA(n_components=152, random_state= 42)
mat_pca = pca.fit_transform(mat_normalized)

6.1.2 K-Means

Elbow Plot

cluster_num_list=range(2,31)
inertia_lis = []
batch_size = 2048

for i in cluster_num_list:
    
    kmeans = MiniBatchKMeans(init='k-means++', n_clusters= i, batch_size=batch_size, n_init=10, 
                             max_no_improvement = 10, verbose=0, random_state=42)
    kmeans.fit(mat_umap_cluster)
        
    inertia_lis.append(kmeans.inertia_)
    print(f"{i} clusters completed.")

d:\miniconda3\envs\pydata-book\Lib\site-packages\joblib\externals\loky\backend\context.py:136: UserWarning:

Could not find the number of physical cores for the following reason:
[WinError 2] 系统找不到指定的文件。
Returning the number of logical cores instead. You can silence this warning by setting LOKY_MAX_CPU_COUNT to the number of cores you want to use.

  File "d:\miniconda3\envs\pydata-book\Lib\site-packages\joblib\externals\loky\backend\context.py", line 257, in _count_physical_cores
    cpu_info = subprocess.run(
               ^^^^^^^^^^^^^^^
  File "d:\miniconda3\envs\pydata-book\Lib\subprocess.py", line 548, in run
    with Popen(*popenargs, **kwargs) as process:
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "d:\miniconda3\envs\pydata-book\Lib\subprocess.py", line 1026, in __init__
    self._execute_child(args, executable, preexec_fn, close_fds,
  File "d:\miniconda3\envs\pydata-book\Lib\subprocess.py", line 1538, in _execute_child
    hp, ht, pid, tid = _winapi.CreateProcess(executable, args,
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

2 clusters completed.
3 clusters completed.
4 clusters completed.
5 clusters completed.
6 clusters completed.
7 clusters completed.
8 clusters completed.
9 clusters completed.
10 clusters completed.
11 clusters completed.
12 clusters completed.
13 clusters completed.
14 clusters completed.
15 clusters completed.
16 clusters completed.
17 clusters completed.
18 clusters completed.
19 clusters completed.
20 clusters completed.
21 clusters completed.
22 clusters completed.
23 clusters completed.
24 clusters completed.
25 clusters completed.
26 clusters completed.
27 clusters completed.
28 clusters completed.
29 clusters completed.
30 clusters completed.

fig, ax = plt.subplots(figsize = (12, 8))
plt.plot(cluster_num_list, inertia_lis)
plt.xlabel('Number of Clusters in K-means')
plt.ylabel('Inertia')
plt.xticks(cluster_num_list, rotation = 45, fontsize = 10)
plt.title('Elbow Plot')
plt.grid(True, alpha=0.3)
plt.show()

No clear elbow displayed, possible optimal number of clusters: 6 and 9.

Calinski-Harabasz Index

calinski_score_list=[]
for k in cluster_num_list:

    kmeans = MiniBatchKMeans(init='k-means++', n_clusters= k, batch_size=batch_size, n_init=10, 
                             max_no_improvement = 10, verbose=0, random_state=42).fit(mat_umap_cluster)
        
    cluster_labels = kmeans.labels_
        
    calinski_score=calinski_harabasz_score(mat_umap_cluster, cluster_labels)
        
    calinski_score_list.append(calinski_score)

    print(f"{k} clusters completed.")

2 clusters completed.
3 clusters completed.
4 clusters completed.
5 clusters completed.
6 clusters completed.
7 clusters completed.
8 clusters completed.
9 clusters completed.
10 clusters completed.
11 clusters completed.
12 clusters completed.
13 clusters completed.
14 clusters completed.
15 clusters completed.
16 clusters completed.
17 clusters completed.
18 clusters completed.
19 clusters completed.
20 clusters completed.
21 clusters completed.
22 clusters completed.
23 clusters completed.
24 clusters completed.
25 clusters completed.
26 clusters completed.
27 clusters completed.
28 clusters completed.
29 clusters completed.
30 clusters completed.

fig, ax = plt.subplots(figsize = (12, 8))
plt.plot(cluster_num_list, calinski_score_list)
plt.xlabel('Number of Clusters in K-means')
plt.xticks(cluster_num_list, rotation = 45, fontsize = 10)
plt.ylabel('Calinski-Harabasz score')
plt.title('Calinski-Harabasz Score Plot')
plt.grid(True, alpha=0.3)
plt.show()

Davies-Bouldin Index.

davies_score_list=[]
for k in cluster_num_list:

    kmeans = MiniBatchKMeans(init='k-means++', n_clusters= k, batch_size=batch_size, n_init=10, 
                             max_no_improvement = 10, verbose=0, random_state=42).fit(mat_umap_cluster)
        
    cluster_labels = kmeans.labels_
        
    davies_score=davies_bouldin_score(mat_umap_cluster, cluster_labels)
        
    davies_score_list.append(davies_score)

    print(f"{k} clusters completed.")

2 clusters completed.
3 clusters completed.
4 clusters completed.
5 clusters completed.
6 clusters completed.
7 clusters completed.
8 clusters completed.
9 clusters completed.
10 clusters completed.
11 clusters completed.
12 clusters completed.
13 clusters completed.
14 clusters completed.
15 clusters completed.
16 clusters completed.
17 clusters completed.
18 clusters completed.
19 clusters completed.
20 clusters completed.
21 clusters completed.
22 clusters completed.
23 clusters completed.
24 clusters completed.
25 clusters completed.
26 clusters completed.
27 clusters completed.
28 clusters completed.
29 clusters completed.
30 clusters completed.

fig, ax = plt.subplots(figsize = (12, 8))
plt.plot(cluster_num_list,davies_score_list)
plt.xlabel('Number of Clusters in K-means')
plt.xticks(cluster_num_list, rotation = 45, fontsize = 10)
plt.ylabel('Davies-Bouldin score')
plt.title('Davies-Bouldin Score Plot')
plt.grid(True, alpha=0.3)
plt.show()

Trustworthiness (based on stratified sample)

umap = UMAP(n_neighbors=45, min_dist = 0.1, n_components=50, init='pca', metric='cosine', n_jobs=-1, verbose=True)
mat_umap_cluster_sample = umap.fit_transform(df_sbert_sample.iloc[:, : -1])

UMAP(angular_rp_forest=True, init='pca', metric='cosine', n_components=50, n_neighbors=45, verbose=True)
Sun Jan 11 00:11:09 2026 Construct fuzzy simplicial set
Sun Jan 11 00:11:09 2026 Finding Nearest Neighbors
Sun Jan 11 00:11:09 2026 Building RP forest with 12 trees
Sun Jan 11 00:11:10 2026 NN descent for 14 iterations
     1  /  14
     2  /  14
     3  /  14
     4  /  14
    Stopping threshold met -- exiting after 4 iterations
Sun Jan 11 00:11:15 2026 Finished Nearest Neighbor Search
Sun Jan 11 00:11:15 2026 Construct embedding

    completed  0  /  200 epochs
    completed  20  /  200 epochs
    completed  40  /  200 epochs
    completed  60  /  200 epochs
    completed  80  /  200 epochs
    completed  100  /  200 epochs
    completed  120  /  200 epochs
    completed  140  /  200 epochs
    completed  160  /  200 epochs
    completed  180  /  200 epochs
Sun Jan 11 00:11:25 2026 Finished embedding

Re-run UMAP

umap = UMAP(n_neighbors=45, min_dist = 0.1, n_components=26, init='pca', metric='cosine', n_jobs=-1, verbose=True)
mat_umap_cluster = umap.fit_transform(df_sbert)

UMAP(angular_rp_forest=True, init='pca', metric='cosine', n_components=26, n_neighbors=45, verbose=True)
Sun Jan 11 00:52:27 2026 Construct fuzzy simplicial set
Sun Jan 11 00:52:28 2026 Finding Nearest Neighbors
Sun Jan 11 00:52:28 2026 Building RP forest with 52 trees
Sun Jan 11 00:53:12 2026 NN descent for 20 iterations
     1  /  20
     2  /  20
     3  /  20
     4  /  20
     5  /  20
    Stopping threshold met -- exiting after 5 iterations
Sun Jan 11 00:56:33 2026 Finished Nearest Neighbor Search
Sun Jan 11 00:56:48 2026 Construct embedding

    completed  0  /  200 epochs
    completed  20  /  200 epochs
    completed  40  /  200 epochs
    completed  60  /  200 epochs
    completed  80  /  200 epochs
    completed  100  /  200 epochs
    completed  120  /  200 epochs
    completed  140  /  200 epochs
    completed  160  /  200 epochs
    completed  180  /  200 epochs
Sun Jan 11 01:05:06 2026 Finished embedding

cluster_num_list=range(2,14)

avg_sil_score_list=[]
for k in cluster_num_list:

    kmeans=KMeans(n_clusters=k, random_state=42).fit(mat_pca_sample)
        
    cluster_labels = kmeans.labels_
        
    avg_sil_score=silhouette_score(mat_pca_sample, cluster_labels)
        
    avg_sil_score_list.append(avg_sil_score)

    print(f"{k} clusters completed.")


plt.plot(cluster_num_list,avg_sil_score_list)
plt.xlabel('Number of Clusters in K-means')
plt.ylabel('Average Silhouette Score')
plt.title('Average Silhouette Score Plot')
plt.grid(True, alpha=0.3)
plt.show()

2 clusters completed.
3 clusters completed.
4 clusters completed.
5 clusters completed.
6 clusters completed.
7 clusters completed.
8 clusters completed.
9 clusters completed.
10 clusters completed.
11 clusters completed.
12 clusters completed.
13 clusters completed.

Notice that the overall average silhouette scores are very low, the silhouette coefficient assumes that clusters are convex and well-separated. However, high-dimensional text embeddings (even dimension-reduced) generated by SBERT often form non-convex, density-based structures where semantic boundaries are fluid, naturally resulting in lower silhouette scores compared to tabular data.

Cluster sorted matrix

cluster_num_list=range(2,14)

for k in cluster_num_list:
    kmeans = KMeans(n_clusters=k, random_state=42).fit(mat_pca)
    
    a = df_sbert_sample.iloc[:, :-1].copy() 
    a['k_mean_label'] = kmeans.labels_
    a_sorted = a.sort_values(by=["k_mean_label"])
    
    print(f'--- K={k} Cluster Distribution ---')
    print(a_sorted['k_mean_label'].value_counts())
    
    data_for_dist = a_sorted.drop(['k_mean_label'], axis=1)
    dist_mat = squareform(pdist(data_for_dist, metric='euclidean'))
    
    fig, ax = plt.subplots(figsize=(10, 8))
    
    c = ax.imshow(dist_mat, cmap='viridis', aspect='auto', origin='lower') 
    
    ax.set_title(f"Distance Matrix Sorted by Cluster (K={k})")
    ax.set_xlabel('Data Points (Sorted)')
    ax.set_ylabel('Data Points (Sorted)')
    fig.colorbar(c, ax=ax, shrink=0.8, label='Euclidean Distance')
    
    plt.show()
    
    plt.close(fig) 

    del dist_mat
    del a_sorted
    del data_for_dist

--- K=2 Cluster Distribution ---
k_mean_label
1    7395
0    2983
Name: count, dtype: int64

--- K=3 Cluster Distribution ---
k_mean_label
1    4147
0    3657
2    2574
Name: count, dtype: int64

--- K=4 Cluster Distribution ---
k_mean_label
1    3818
3    3477
0    1990
2    1093
Name: count, dtype: int64

--- K=5 Cluster Distribution ---
k_mean_label
1    3640
3    3286
0    1942
2    1087
4     423
Name: count, dtype: int64

--- K=6 Cluster Distribution ---
k_mean_label
1    3090
0    2411
5    1801
3    1640
2    1036
4     400
Name: count, dtype: int64

--- K=7 Cluster Distribution ---
k_mean_label
6    2847
0    2294
5    1734
3    1086
1    1025
2     995
4     397
Name: count, dtype: int64

--- K=8 Cluster Distribution ---
k_mean_label
5    2069
7    1745
0    1716
4    1622
3    1308
2     986
1     537
6     395
Name: count, dtype: int64

--- K=9 Cluster Distribution ---
k_mean_label
7    2507
5    2168
4    1297
1     999
3     995
2     945
0     549
8     517
6     401
Name: count, dtype: int64

--- K=10 Cluster Distribution ---
k_mean_label
5    1701
4    1534
9    1497
7    1309
1    1023
2     967
3     864
0     662
8     441
6     380
Name: count, dtype: int64

--- K=11 Cluster Distribution ---
k_mean_label
5     1708
9     1442
7     1315
0     1195
3      936
2      915
10     859
1      655
4      536
8      439
6      378
Name: count, dtype: int64

--- K=12 Cluster Distribution ---
k_mean_label
5     1594
9     1441
7     1409
0     1186
2      913
10     858
1      648
3      561
4      528
8      439
11     423
6      378
Name: count, dtype: int64

--- K=13 Cluster Distribution ---
k_mean_label
5     1264
12    1166
0     1156
7     1086
3     1037
2      904
10     781
9      623
1      580
4      525
8      444
11     432
6      380
Name: count, dtype: int64

As the number of clusters \(k\) increases, a distinct block-diagonal structure emerges. The diagonal regions gradually become darker, indicating improved intra-cluster cohesion, while the off-diagonal areas become lighter, signaling better inter-cluster separation. However, if we keep increase \(k\), we observe a degradation in separation, where clusters begin to exhibit high similarity to one another again.

Above three clustering metrics presented, we choose \(k\) = 9.

6.3 Re-run K-Means

batch_size = 2048
kmeans = MiniBatchKMeans(init='k-means++', n_clusters= 5, batch_size=batch_size, n_init=10, 
                         max_no_improvement = 10, verbose=0, random_state=42).fit(mat_umap_cluster)

Color code by predicted k-means labels.

df_umap['k_mean_label'] = kmeans.labels_


fig, ax = plt.subplots()
sns.scatterplot(x=df_umap.iloc[:,0], y=df_umap.iloc[:,1], hue=df_umap['k_mean_label'])
ax.set_title(f"Perplexity Value {40} and Random State {42}")
ax.set_xlabel("x_projected")
ax.set_ylabel("y_projected")
plt.show()

We see the clustering result is not perfect, as some isolated clusters contain multiple k-means cluster labels. However, noticing the algorithm is splitting a larger cluster into multiple closed sub-clusters, which is indeed indicating that the algorithm is capturing granular relationships. So the clustering result is somewhat reasonable.

fig, ax = plt.subplots()
sns.scatterplot(x=df_umap.iloc[:,0], y=df_umap.iloc[:,1], color = 'lightgray', alpha = 0.4)

mask = df_umap['k_mean_label'] == 2
sns.scatterplot(x=df_umap[mask].iloc[:,0], y=df_umap[mask].iloc[:,1], color = 'red', alpha = 0.9)

ax.set_title(f"Perplexity Value {40} and Random State {42}")
ax.set_xlabel("x_projected")
ax.set_ylabel("y_projected")
plt.show()

6.4 PCA + DBSCAN

The input matrix is also using the same PCA-dimension-reduced matrix.

Plot K-distance plot for k = 2 ~ 10:

#Find the distance of the kth nearest neighbor
for k in range(2,11):
    nbrs = NearestNeighbors(n_neighbors=k, algorithm='ball_tree').fit(mat_pca)
    distances, indices = nbrs.kneighbors(mat_pca)

    plt.plot(np.sort(distances.T[k-1]))
    plt.xlabel('Points Sorted by Distance to k=%s Nearest Neighbor'%k)
    plt.ylabel('k=%s Nearest Neighbor Distance'%k)
    plt.title('Noisy Dataset')
    plt.show()

Based on the k-distance plot with \(k=10\), we identified the most distinct elbow point around the distance value of 0.78. The region below this threshold represents the dense clusters, while the sharp rise beyond 0.8 indicates outliers. Therefore, we selected eps=0.78 and min_pts=10 as our baseline parameters for DBSCAN to balance cluster granularity and noise reduction.

We set the patameter range to \(\epsilon \in\) [0.5, 0.95] with a step of 0.05, min_pts \(\in\) [5, 10, 15, 20] to do a small range best parameter search.

X_input = mat_pca 

eps_range = np.arange(0.5, 0.95, 0.05) 

min_samples_list = [5, 10, 15, 20] 

for minpts in min_samples_list:
    silhouette_scores = []
    num_clusters_list = []
    num_noise_points = []
    valid_eps = [] 

    print(f"--- Processing min_samples = {minpts} ---")

    for eps in eps_range:
        dst = DBSCAN(eps=eps, min_samples=minpts, metric='euclidean', n_jobs=-1)
        labels = dst.fit_predict(X_input)

        n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
        n_noise = list(labels).count(-1)
        
        num_clusters_list.append(n_clusters)
        num_noise_points.append(n_noise)
        valid_eps.append(eps)
        
        if n_clusters > 1 and n_noise < len(X_input):
            non_noise_mask = labels != -1
            
            score = silhouette_score(X_input[non_noise_mask], labels[non_noise_mask])
            silhouette_scores.append(score)
        else:
            silhouette_scores.append(0)
        
    fig, ax = plt.subplots(1, 3, figsize=(18, 5), constrained_layout=True)
    fig.suptitle(f'DBSCAN Parameter Search (min_samples={minpts})', fontsize=16)

    ax[0].plot(valid_eps, silhouette_scores, marker='o', color='tab:blue')
    ax[0].set_title('Avg Silhouette Score')
    ax[0].set_xlabel('Epsilon')
    ax[0].set_ylabel('Score')
    ax[0].grid(True, linestyle='--', alpha=0.6)

    ax[1].plot(valid_eps, num_noise_points, marker='o', color='tab:red')
    ax[1].set_title('Number of Noise Points')
    ax[1].set_xlabel('Epsilon')
    ax[1].set_ylabel('Count')
    ax[1].grid(True, linestyle='--', alpha=0.6)

    ax[2].plot(valid_eps, num_clusters_list, marker='o', color='tab:green')
    ax[2].set_title('Number of Clusters')
    ax[2].set_xlabel('Epsilon')
    ax[2].set_ylabel('Count')
    ax[2].grid(True, linestyle='--', alpha=0.6)

    plt.show()

--- Processing min_samples = 5 ---

--- Processing min_samples = 10 ---

--- Processing min_samples = 15 ---

--- Processing min_samples = 20 ---

While the grid search showed highest Silhouette scores at \(\epsilon \ge 0.70\), these settings resulted in too little noise detected, and meaningful clusters are merged into 1-2 larger groups. We identified the optimal operating point at \(\epsilon =0.65\) with min_samples=5.

6.5 Re-run DBSCAN + Clustering Analysis

BEST_EPS = 0.65        
BEST_MIN_SAMPLES = 5   
X_INPUT = mat_pca

db = DBSCAN(eps=BEST_EPS, min_samples=BEST_MIN_SAMPLES, metric='euclidean', n_jobs=-1)
labels = db.fit_predict(X_INPUT)

a = df_sbert_sample.iloc[:, :-1].copy() 
a['db_label'] = labels
a_sorted = a.sort_values(by=["db_label"])

print(a_sorted['db_label'].value_counts())

data_for_dist = a_sorted.drop(['db_label'], axis=1)
dist_mat = squareform(pdist(data_for_dist, metric='euclidean'))


fig, ax = plt.subplots(figsize=(12, 10))

im = ax.imshow(dist_mat, cmap='viridis', aspect='auto', origin='lower')

cbar = plt.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('Euclidean Distance')


ax.set_title(f"DBSCAN Sorted Distance Matrix\n(eps={BEST_EPS}, min_samples={BEST_MIN_SAMPLES}")
ax.set_xlabel('Data Points (Sorted by Cluster Label)')
ax.set_ylabel('Data Points (Sorted by Cluster Label)')

plt.tight_layout()
plt.show()

del dist_mat

db_label
 0    7547
-1    2777
 4      26
 1      10
 2       7
 3       7
 5       4
Name: count, dtype: int64

df_tsne_sample['db_label'] = labels


fig, ax = plt.subplots()
sns.scatterplot(x=df_tsne_sample.iloc[:,0], y=df_tsne_sample.iloc[:,1], hue=df_tsne_sample['db_label'])
ax.set_title(f"Perplexity Value {40} and Random State {42}")
ax.set_xlabel("x_projected")
ax.set_ylabel("y_projected")
plt.show()

However, by visually checking the post-cluster sorted matrix and the color-coded tsne plot by predicted label, the performance of DBSCAN was clearly worse that of K-Means. The clusters appeared less distinct and contained misclasscified noise. Consequently, we choose K-Means with \(k=9\) as our optimal clustering solution.

7. Interpretate Result

kmeans = KMeans(n_clusters= 9, random_state=42)
kmeans.fit(mat_pca)

KMeans(n_clusters=9, random_state=42)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

df_sample = review_df_all_deduplicate.loc[sample_idx]
df_sample['predicted_label'] = kmeans.labels_

Plot all the numercial varaibles (except massive na valued variables) across clusters.

num_col = ['rating', 'total_feedback_count', 'total_neg_feedback_count', 'total_pos_feedback_count']

for num in num_col :
    plt.figure(figsize=(6,4))
    sns.boxplot(data=df_sample, x='predicted_label', y='rating')  
    plt.xlabel('Predicted cluster')
    plt.ylabel(f'{num}')
    plt.title('num_col by cluster')
    plt.tight_layout()
    plt.show()

print(df_sample.groupby('predicted_label').agg({'rating': 'mean', 'review_text': 'count'    
}).sort_values(by='rating'))

                   rating  review_text
predicted_label                       
1                3.994995          999
3                4.013065          995
7                4.183486         2507
6                4.359102          401
8                4.361702          517
4                4.387047         1297
5                4.399446         2168
2                4.400000          945
0                4.520947          549

The table summarizes the statistical profile of the 9 identified clusters, sorted by average rating. We observe a narrow rating distribution ranging from 3.99 to 4.52, indicating an overall positive sentiment across the sampled dataset. Cluster 7 is the dominant group by volume, while Cluster 1 and Cluster 0 represent the lowest and highest satisfaction extremes, respectively.

sorted_clusters = sorted(df_sample['predicted_label'].unique())


print("-" * 60)

for cluster_id in sorted_clusters:

    cluster_data = df_sample[df_sample['predicted_label'] == cluster_id]
    
    combined_text = " ".join(cluster_data['review_text_cleaned'].astype(str).tolist())

    ngrams_list = generate_ngrams(combined_text, n_gram=2)
    
    top_keywords = Counter(ngrams_list).most_common(10)
    
    print(f"Cluster {cluster_id} (Size: {len(cluster_data)})")
    print(f"Top Keywords: {top_keywords}")
    print("-" * 60)

------------------------------------------------------------
Cluster 0 (Size: 549)
Top Keywords: [('dermalogica sampling', 91), ('gifted dermalogica', 78), ('dead skin', 73), ('sensitive skin', 72), ('daily microfoliant', 59), ('skin feeling', 52), ('complimentary dermalogica', 52), ('leaves skin', 49), ('skin feels', 49), ('dry skin', 44)]
------------------------------------------------------------
Cluster 1 (Size: 999)
Top Keywords: [('full size', 68), ('love product', 51), ('received product', 40), ('received sample', 35), ('sample size', 33), ('skin feel', 31), ('goes long', 29), ('long way', 29), ('highly recommend', 27), ('product free', 27)]
------------------------------------------------------------
Cluster 2 (Size: 945)
Top Keywords: [('makeup remover', 156), ('eye makeup', 129), ('remove makeup', 94), ('removes makeup', 78), ('waterproof mascara', 77), ('cleansing balm', 76), ('leaves skin', 69), ('skin feeling', 62), ('removing makeup', 60), ('love product', 58)]
------------------------------------------------------------
Cluster 3 (Size: 995)
Top Keywords: [('received product', 69), ('leaves skin', 60), ('sensitive skin', 59), ('skin feel', 56), ('long way', 55), ('goes long', 53), ('love product', 48), ('dry skin', 47), ('makes skin', 47), ('skin feeling', 46)]
------------------------------------------------------------
Cluster 4 (Size: 1297)
Top Keywords: [('dry skin', 137), ('sensitive skin', 126), ('leaves skin', 115), ('face wash', 114), ('skin feeling', 102), ('cleansing balm', 100), ('love cleanser', 97), ('skin feels', 91), ('oily skin', 81), ('acne prone', 64)]
------------------------------------------------------------
Cluster 5 (Size: 2168)
Top Keywords: [('dry skin', 374), ('skin feels', 171), ('sensitive skin', 166), ('skin feel', 151), ('oily skin', 143), ('makes skin', 136), ('long way', 122), ('goes long', 118), ('received product', 112), ('skin feeling', 109)]
------------------------------------------------------------
Cluster 6 (Size: 401)
Top Keywords: [('dry skin', 48), ('jet lag', 47), ('summer fridays', 42), ('love mask', 42), ('lag mask', 39), ('sensitive skin', 36), ('overnight mask', 27), ('face mask', 23), ('makes skin', 19), ('next morning', 18)]
------------------------------------------------------------
Cluster 7 (Size: 2507)
Top Keywords: [('sensitive skin', 269), ('acne prone', 162), ('love product', 160), ('using product', 152), ('prone skin', 132), ('dry skin', 119), ('difference skin', 117), ('skin tone', 116), ('skin feels', 113), ('made skin', 111)]
------------------------------------------------------------
Cluster 8 (Size: 517)
Top Keywords: [('sun protection', 154), ('protection factor', 136), ('white cast', 66), ('love sunscreen', 39), ('best sunscreen', 32), ('oily skin', 24), ('sunscreen used', 24), ('leave white', 23), ('doesnt leave', 23), ('unseen sunscreen', 23)]
------------------------------------------------------------

The largest portion of the dataset is driven by efficacy regarding specific skin types. Cluster 7 (the largest group) and Cluster 4 focus heavily on problem-solving for “acne prone” and “sensitive skin,” indicating that efficacy is the primary driver of engagement. Meanwhile, Cluster 5 captures the positive experiences of “dry skin” users seeking hydration.

The model isolated distinct product lines. Cluster 8 identifies “sun protection” with a specific negative sentiment regarding “white cast,” highlighting a critical product defect. Cluster 6 is unique as it is dominated by a single viral product (“Summer Fridays Jet Lag Mask”), demonstrating how specific ‘hero products’ can form their own semantic clusters.

While cluster 0, 1, and 3 are the non-product-related clusters. Cluster 0 is characterized by terms like “gifted” and “dermalogica sampling,” representing incentivized reviews that may introduce positive bias. Conversely, Cluster 1 focuses on “sample size” complaints.

Visualization of top 50 bigrams for each cluster.

n_clusters = 9
n_cols = 3
n_rows = math.ceil(n_clusters / n_cols)

fig, axes = plt.subplots(n_rows, n_cols, figsize=(20, 15))
axes = axes.flatten() 

sorted_clusters = sorted(df_sample['predicted_label'].unique())


for i, cluster_id in enumerate(sorted_clusters):
    ax = axes[i]
    
    cluster_data = df_sample[df_sample['predicted_label'] == cluster_id]
    
    combined_text = " ".join(cluster_data['review_text_cleaned'].astype(str).tolist())
    
    ngrams_list = generate_ngrams(combined_text, n_gram=2)
    
    ngram_counts = dict(Counter(ngrams_list).most_common(50))
    
    if len(ngram_counts) > 0:
        wc = WordCloud(
            width=800, 
            height=400, 
            background_color='white', 
            colormap='viridis', 
            max_font_size=100
        ).generate_from_frequencies(ngram_counts)
        
        ax.imshow(wc, interpolation="bilinear")
        ax.set_title(f"Cluster {cluster_id} Theme\n(Size: {len(cluster_data)})", fontsize=14)
    else:
        ax.text(0.5, 0.5, "Not enough data", ha='center', va='center')
        
    ax.axis('off') 

for j in range(i + 1, len(axes)):
    axes[j].axis('off')

plt.tight_layout()
plt.show()

8. Conclusion

To conclude, by replacing static star ratings with an unsupervised learning result, we transformed a sample unstructured customer feedback into 9 distinct semantic clusters. We experimented two types of algorithms: K-Means and DBSCAN, and by comparison of the 2 algorithms, we found the K-Means is a better fit to our case. Then we interpreted our clustering results in terms of their business value, and classified the 9 clusterings into 3 distinct categories: skin-related topic, product-related topic, and non-product-related noise topic.

Next Step

In this project, the whole clustering building process is only using around 10,000 sample reviews, the clustering result may not represent the whole reviews. Additionally, we identify the limitation of methods we used in this project, such as the computation limitation of the t-sne, and not the most ideal cluster assumptions for both algorithms, a more robust approach is needed.